1 2 3 #include "aij.h" 4 #include "inline/spops.h" 5 /* 6 Factorization code for AIJ format. 7 */ 8 9 int MatiAIJLUFactorSymbolic(Mat mat,IS isrow,IS iscol,Mat *fact) 10 { 11 Matiaij *aij = (Matiaij *) mat->data, *aijnew; 12 IS isicol; 13 int *r,*ic, ierr, i, j, n = aij->m, *ai = aij->i, *aj = aij->j; 14 int prow, *ainew,*ajnew, jmax,*fill, *ajtmp, nz , *ii; 15 int *idnew, idx, pivot_row,row,m,fm, nnz, nzi,len; 16 17 if (n != aij->n) SETERR(1,"Mat must be square"); 18 if (!isrow) {SETERR(1,"Must have row permutation");} 19 if (!iscol) {SETERR(1,"Must have column permutation");} 20 21 if (ierr = ISInvertPermutation(iscol,&isicol)) SETERR(ierr,0); 22 ISGetIndices(isrow,&r); ISGetIndices(isicol,&ic); 23 24 /* get new row pointers */ 25 ainew = (int *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(ainew); 26 ainew[0] = 1; 27 /* don't know how many column pointers are needed so estimate */ 28 jmax = 2*ai[n]; 29 ajnew = (int *) MALLOC( (jmax)*sizeof(int) ); CHKPTR(ajnew); 30 /* fill is a linked list of nonzeros in active row */ 31 fill = (int *) MALLOC( (n+1)*sizeof(int)); CHKPTR(fill); 32 /* idnew is location of diagonal in factor */ 33 idnew = (int *) MALLOC( (n+1)*sizeof(int)); CHKPTR(idnew); 34 idnew[0] = 1; 35 36 for ( i=0; i<n; i++ ) { 37 /* first copy previous fill into linked list */ 38 nnz = nz = ai[r[i]+1] - ai[r[i]]; 39 ajtmp = aj + ai[r[i]] - 1; 40 fill[n] = n; 41 while (nz--) { 42 fm = n; 43 idx = ic[*ajtmp++ - 1]; 44 do { 45 m = fm; 46 fm = fill[m]; 47 } while (fm < idx); 48 fill[m] = idx; 49 fill[idx] = fm; 50 } 51 row = fill[n]; 52 while ( row < i ) { 53 ajtmp = ajnew + idnew[row] - 1; 54 nz = ainew[row+1] - idnew[row]; 55 fm = row; 56 while (nz--) { 57 fm = n; 58 idx = *ajtmp++ - 1; 59 do { 60 m = fm; 61 fm = fill[m]; 62 } while (fm < idx); 63 if (fm != idx) { 64 fill[m] = idx; 65 fill[idx] = fm; 66 fm = idx; 67 nnz++; 68 } 69 } 70 row = fill[row]; 71 } 72 /* copy new filled row into permanent storage */ 73 ainew[i+1] = ainew[i] + nnz; 74 if (ainew[i+1] > jmax+1) { 75 /* allocate a longer ajnew */ 76 jmax += nnz*(n-i); 77 ajtmp = (int *) MALLOC( jmax*sizeof(int) );CHKPTR(ajtmp); 78 MEMCPY(ajtmp,ajnew,(ainew[i]-1)*sizeof(int)); 79 FREE(ajnew); 80 ajnew = ajtmp; 81 } 82 ajtmp = ajnew + ainew[i] - 1; 83 fm = fill[n]; 84 nzi = 0; 85 while (nnz--) { 86 if (fm < i) nzi++; 87 *ajtmp++ = fm + 1; 88 fm = fill[fm]; 89 } 90 idnew[i] = ainew[i] + nzi; 91 } 92 93 ISDestroy(isicol); FREE(fill); 94 95 /* put together the new matrix */ 96 ierr = MatCreateSequentialAIJ(n, n, 0, fact); CHKERR(ierr); 97 aijnew = (Matiaij *) (*fact)->data; 98 FREE(aijnew->imax); 99 aijnew->singlemalloc = 0; 100 len = (ainew[n] - 1)*sizeof(double); 101 aijnew->a = (Scalar *) MALLOC( len ); CHKPTR(aijnew->a); 102 aijnew->j = ajnew; 103 aijnew->i = ainew; 104 aijnew->diag = idnew; 105 (*fact)->row = isrow; 106 (*fact)->col = iscol; 107 (*fact)->factor = FACTOR_LU; 108 return 0; 109 } 110 111 int MatiAIJLUFactorNumeric(Mat mat,Mat fact) 112 { 113 Matiaij *aij = (Matiaij *) mat->data, *aijnew = (Matiaij *)fact->data; 114 IS iscol = fact->col, isrow = fact->row, isicol; 115 int *r,*ic, ierr, i, j, n = aij->m, *ai = aijnew->i, *aj = aijnew->j; 116 int prow, *ainew,*ajnew, jmax,*fill, *ajtmpold, *ajtmp, nz , *ii; 117 int *idnew, idx, pivot_row,row,*pj, m,fm, nnz, nzi,len; 118 Scalar *rtmp,*vnew,*v, *pv, *pc, multiplier; 119 120 if (ierr = ISInvertPermutation(iscol,&isicol)) SETERR(ierr,0); 121 ierr = ISGetIndices(isrow,&r); CHKERR(ierr); 122 ierr = ISGetIndices(isicol,&ic); CHKERR(ierr); 123 rtmp = (Scalar *) MALLOC( (n+1)*sizeof(Scalar) ); CHKPTR(rtmp); 124 125 for ( i=0; i<n; i++ ) { 126 nz = ai[i+1] - ai[i]; 127 ajtmp = aj + ai[i] - 1; 128 for ( j=0; j<nz; j++ ) rtmp[ajtmp[j]-1] = 0.0; 129 130 /* load in initial (unfactored row) */ 131 nz = aij->i[r[i]+1] - aij->i[r[i]]; 132 ajtmpold = aij->j + aij->i[r[i]] - 1; 133 v = aij->a + aij->i[r[i]] - 1; 134 for ( j=0; j<nz; j++ ) rtmp[ic[ajtmpold[j]-1]] = v[j]; 135 136 row = *ajtmp++ - 1; 137 while (row < i) { 138 pc = rtmp + row; 139 if (*pc != 0.0) { 140 nz = aijnew->diag[row] - ai[row]; 141 pv = aijnew->a + aijnew->diag[row] - 1; 142 pj = aijnew->j + aijnew->diag[row]; 143 multiplier = *pc * *pv++; 144 *pc = multiplier; 145 nz = ai[row+1] - ai[row] - 1 - nz; 146 while (nz-->0) rtmp[*pj++ - 1] -= multiplier* *pv++; 147 } 148 row = *ajtmp++ - 1; 149 } 150 /* finished row so stick it into aijnew->a */ 151 pv = aijnew->a + ai[i] - 1; 152 pj = aijnew->j + ai[i] - 1; 153 nz = ai[i+1] - ai[i]; 154 rtmp[i] = 1.0/rtmp[i]; 155 for ( j=0; j<nz; j++ ) {pv[j] = rtmp[pj[j]-1];} 156 } 157 FREE(rtmp); 158 ierr = ISDestroy(isicol); CHKERR(ierr); 159 fact->factor = FACTOR_LU; 160 161 return 0; 162 } 163 int MatiAIJSolve(Mat mat,Vec bb, Vec xx) 164 { 165 Matiaij *aij = (Matiaij *) mat->data; 166 IS iscol = mat->col, isrow = mat->row; 167 int *r,*c, ierr, i, j, n = aij->m, *vi, *ai = aij->i, *aj = aij->j; 168 int nz; 169 Scalar *x,*b,*tmp, *aa = aij->a, sum, *v; 170 171 if (ierr = VecGetArray(bb,&b)) SETERR(ierr,0); 172 if (ierr = VecGetArray(xx,&x)) SETERR(ierr,0); 173 tmp = (Scalar *) MALLOC(n*sizeof(Scalar)); CHKPTR(tmp); 174 175 if (ierr = ISGetIndices(isrow,&r)) SETERR(ierr,0); 176 if (ierr = ISGetIndices(iscol,&c)) SETERR(ierr,0); c = c + (n-1); 177 178 /* forward solve the lower triangular */ 179 tmp[0] = b[*r++]; 180 for ( i=1; i<n; i++ ) { 181 v = aa + ai[i] - 1; 182 vi = aj + ai[i] - 1; 183 nz = aij->diag[i] - ai[i]; 184 sum = b[*r++]; 185 while (nz--) sum -= *v++ * tmp[*vi++ - 1]; 186 tmp[i] = sum; 187 } 188 189 /* backward solve the upper triangular */ 190 for ( i=n-1; i>=0; i-- ){ 191 v = aa + aij->diag[i]; 192 vi = aj + aij->diag[i]; 193 nz = ai[i+1] - aij->diag[i] - 1; 194 sum = tmp[i]; 195 while (nz--) sum -= *v++ * tmp[*vi++ - 1]; 196 x[*c--] = tmp[i] = sum*aa[aij->diag[i]-1]; 197 } 198 199 FREE(tmp); 200 return 0; 201 } 202