1 #include <omp.h>
    2 
    3 #include "matrix.h"
    4 #include <math.h>
    5 
    6 double *TriangularSolve(L, b)
    7 SMatrix L;
    8 double *b;
    9 {
   10         int i, j;
   11         int *theRow;
   12         double *y, *x;
   13 
   14         x = NewVector(L.n);
   15         y = NewVector(L.n);
   16 
   17         /* forward solve */
   18 
   19         for (j=0; j<L.n; j++) {
   20                 theRow = L.row + (L.startrow[j]-L.firstnz[j]);
   21                 y[j] = b[j]/L.nz[L.firstnz[j]];
   22                 for (i=L.firstnz[j]+1; i<L.lastnz[j]; i++)
   23                         b[theRow[i]] -= L.nz[i]*y[j];
   24                 }
   25 
   26         /* back solve */
   27 
   28    #pragma omp parallel for private (i, theRow) schedule(dynamic)
   29         for (j=L.n-1; j>=0; j--) {
   30                 theRow = L.row + (L.startrow[j]-L.firstnz[j]);
   31                 for (i=L.firstnz[j]+1; i<L.lastnz[j]; i++)
   32                         y[j] -= L.nz[i]*x[theRow[i]];
   33                 x[j] = y[j]/L.nz[L.firstnz[j]];
   34                 }
   35 
   36         return(x);
   37 }
   38 
   39 double ComputeNorm(x, n)
   40 double *x;
   41 {
   42         double tmp = 0.0;
   43         int i;
   44 
   45         for (i=0; i<n; i++)
   46                 if (fabs(x[i]-1.0) > tmp)
   47                         tmp = fabs(x[i]-1.0);
   48 
   49         return(tmp);
   50 }
   51 
   52 double *CreateVector(M)
   53 SMatrix M;
   54 {
   55         int i, j;
   56         double tmp, *b;
   57 
   58         b = NewVector(M.n);
   59 
   60         for (j=0; j<M.n; j++)
   61                 for (i=M.firstnz[j]; i<M.firstnz[j+1]; i++) {
   62                         b[M.row[i]] += M.nz[i];
   63 			if (M.row[i] != j)
   64 			  b[j] += M.nz[i];
   65 		      }
   66 
   67         return(b);
   68 }