mthrdslv.c File Reference

#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include "mthrdslv.h"

Include dependency graph for mthrdslv.c:

Go to the source code of this file.

Macros
#define	BIGINT   (1 << 30)
#define	MINPIV   1.0e-8

Typedefs
typedef long int	integer
typedef double	doublereal

Functions
int	naivfct (RMAT a, integer neq, integer *nzr, IMAT ri, integer *nzc, IMAT ci, NZMAT nz, integer *piv, doublereal minpiv)
int	naivslv (RMAT a, integer neq, integer *nzc, IMAT ci, doublereal *rhs, doublereal *sol, integer *piv)

Macro Definition Documentation

#define BIGINT   (1 << 30)

Definition at line 71 of file mthrdslv.c.

#define MINPIV   1.0e-8

Definition at line 73 of file mthrdslv.c.

Referenced by naivfct().

Typedef Documentation

typedef double doublereal

Definition at line 63 of file mthrdslv.c.

typedef long int integer

Definition at line 62 of file mthrdslv.c.

Function Documentation

int naivfct	(	RMAT	a,
		integer	neq,
		integer *	nzr,
		IMAT	ri,
		integer *	nzc,
		IMAT	ci,
		NZMAT	nz,
		integer *	piv,
		doublereal	minpiv
	)

Definition at line 76 of file mthrdslv.c.

References grad::fabs(), MINPIV, NAIVE_ENOPIV, NAIVE_ENULCOL, NAIVE_ERANGE, and NAIVE_MAX.

Referenced by NaiveSolver::Factor().

{
        char todo[neq];
        integer i, j, k, pvr, pvc, nr, nc, r;
        integer *pri, *pci;
        char *prik;
        doublereal den, mul, mulpiv, fari;
#if 0
        doublereal fapvr;
#endif
        doublereal *par, *papvr;

        if (neq <= 0 || (unsigned long)neq > NAIVE_MAX) {
                return NAIVE_ERANGE;
        }

        if (!minpiv) {
                minpiv = MINPIV;
        }
        for (pvr = 0; pvr < neq; pvr++) {
                todo[pvr] = 1;
        }
        for (i = 0; i < neq; i++) {
                if (!nzr[i]) { return NAIVE_ENULCOL + i; }
                nc = neq + 1;   
                nr = nzr[i];
                mul = 0.0;
                pri = ri[i];
                pvr = pri[0];
                mulpiv = 0.;
#if 0
                fapvr = 0.;
                for (k = 0; k < nr; k++) {
                        r = pri[k];
                        if (todo[r]) {
                                fari = fabs(a[r][i]);
                                if (fari > mul) {
                                        mulpiv = fari*minpiv;
                                        if (nzc[r] <= nc  || mulpiv > fapvr) {
                                                nc = nzc[pvr = r];
                                                fapvr = fari;
                                        }
                                        mul = fari;
                                } else if (nzc[r] < nc && fari > mulpiv) {
                                        nc = nzc[pvr = r];
                                }
                        }
                }
#endif
                for (k = 0; k < nr; k++) {
                        r = pri[k];
                        if (todo[r]) {
                                fari = fabs(a[r][i]);
                                if (fari > mul) {
                                        mul = fari;
                                }
                        }
                }
                mulpiv = mul*minpiv;
                for (k = 0; k < nr; k++) {
                        r = pri[k];
                        if (todo[r]) {
                                fari = fabs(a[r][i]);
                                if (fari >= mulpiv && nzc[r] < nc) {
                                        nc = nzc[pvr = r];
                                }
                        }
                }
                if (nc == neq + 1) { return NAIVE_ENOPIV + i; }
                if (mulpiv == 0.)  { return NAIVE_ENOPIV + i; }

                piv[i] = pvr;
                todo[pvr] = 0;
                papvr = a[pvr];
                den = papvr[i] = 1.0/papvr[i];

                for (k = 0; k < nr; k++) {
                        if (!todo[r = pri[k]]) { continue; }
                        par = a[r];
                        mul = par[i] = par[i]*den;
                        prik = nz[r];
                        pci = ci[pvr];
                        for (j = 0; j < nc; j++) {
                                if ((pvc = pci[j]) <= i) { continue; }
                                if (prik[pvc]) {
                                        par[pvc] -= mul*papvr[pvc];
                                } else {
                                        par[pvc] = -mul*papvr[pvc];
                                        prik[pvc] = 1;
                                        ri[pvc][nzr[pvc]++] = r;
                                        ci[r][nzc[r]++] = pvc;
                                }
                        }
                }
        }
        return 0;
}

Here is the call graph for this function:

int naivslv	(	RMAT	a,
		integer	neq,
		integer *	nzc,
		IMAT	ci,
		doublereal *	rhs,
		doublereal *	sol,
		integer *	piv
	)

Definition at line 190 of file mthrdslv.c.

References c, NAIVE_ERANGE, and NAIVE_MAX.

Referenced by NaiveSolver::Solve().

{
        doublereal fwd[neq];

        integer i, k, nc, r, c;
        integer *pci;
        doublereal s;
        doublereal *par;

        if (neq <= 0 || (unsigned long)neq > NAIVE_MAX) {
                return NAIVE_ERANGE;
        }

        fwd[0] = rhs[piv[0]];
        for (i = 1; i < neq; i++) {
                nc = nzc[r = piv[i]];
                s = rhs[r];
                par = a[r];
                pci = ci[r];
                for (k = 0; k < nc; k++) {
                        if ((c = pci[k]) < i) {
                                s -= par[c]*fwd[c];
                        }
                }
                fwd[i] = s;
        }

        r = piv[--neq];
        sol[neq] = fwd[neq]*a[r][neq];
        for (i = neq - 1; i >= 0; i--) {
                nc = nzc[r = piv[i]];
                s = fwd[i];
                par = a[r];
                pci = ci[r];
                for (k = 0; k < nc; k++) {
                        if ((c = pci[k]) > i) {
                                s -= par[c]*sol[c];
                        }
                }
                sol[i] = s*par[i];
        }

        return 0;
}