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 }