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, 0, fact); CHKERR(ierr); 97 aijnew = (Matiaij *) (*fact)->data; 98 FREE(aijnew->imax); 99 aijnew->singlemalloc = 0; 100 len = (ainew[n] - 1)*sizeof(Scalar); 101 /* the next line frees the default space generated by the Create() */ 102 FREE(aijnew->a); FREE(aijnew->ilen); 103 aijnew->a = (Scalar *) MALLOC( len ); CHKPTR(aijnew->a); 104 aijnew->j = ajnew; 105 aijnew->i = ainew; 106 aijnew->diag = idnew; 107 aijnew->ilen = 0; 108 aijnew->imax = 0; 109 (*fact)->row = isrow; 110 (*fact)->col = iscol; 111 (*fact)->factor = FACTOR_LU; 112 return 0; 113 } 114 115 int MatiAIJLUFactorNumeric(Mat mat,Mat *infact) 116 { 117 Mat fact = *infact; 118 Matiaij *aij = (Matiaij *) mat->data, *aijnew = (Matiaij *)fact->data; 119 IS iscol = fact->col, isrow = fact->row, isicol; 120 int *r,*ic, ierr, i, j, n = aij->m, *ai = aijnew->i, *aj = aijnew->j; 121 int prow, *ainew,*ajnew, jmax,*fill, *ajtmpold, *ajtmp, nz , *ii; 122 int *idnew, idx, pivot_row,row,*pj, m,fm, nnz, nzi,len; 123 Scalar *rtmp,*vnew,*v, *pv, *pc, multiplier; 124 125 if (ierr = ISInvertPermutation(iscol,&isicol)) SETERR(ierr,0); 126 ierr = ISGetIndices(isrow,&r); CHKERR(ierr); 127 ierr = ISGetIndices(isicol,&ic); CHKERR(ierr); 128 rtmp = (Scalar *) MALLOC( (n+1)*sizeof(Scalar) ); CHKPTR(rtmp); 129 130 for ( i=0; i<n; i++ ) { 131 nz = ai[i+1] - ai[i]; 132 ajtmp = aj + ai[i] - 1; 133 for ( j=0; j<nz; j++ ) rtmp[ajtmp[j]-1] = 0.0; 134 135 /* load in initial (unfactored row) */ 136 nz = aij->i[r[i]+1] - aij->i[r[i]]; 137 ajtmpold = aij->j + aij->i[r[i]] - 1; 138 v = aij->a + aij->i[r[i]] - 1; 139 for ( j=0; j<nz; j++ ) rtmp[ic[ajtmpold[j]-1]] = v[j]; 140 141 row = *ajtmp++ - 1; 142 while (row < i) { 143 pc = rtmp + row; 144 if (*pc != 0.0) { 145 nz = aijnew->diag[row] - ai[row]; 146 pv = aijnew->a + aijnew->diag[row] - 1; 147 pj = aijnew->j + aijnew->diag[row]; 148 multiplier = *pc * *pv++; 149 *pc = multiplier; 150 nz = ai[row+1] - ai[row] - 1 - nz; 151 while (nz-->0) rtmp[*pj++ - 1] -= multiplier* *pv++; 152 } 153 row = *ajtmp++ - 1; 154 } 155 /* finished row so stick it into aijnew->a */ 156 pv = aijnew->a + ai[i] - 1; 157 pj = aijnew->j + ai[i] - 1; 158 nz = ai[i+1] - ai[i]; 159 rtmp[i] = 1.0/rtmp[i]; 160 for ( j=0; j<nz; j++ ) {pv[j] = rtmp[pj[j]-1];} 161 } 162 FREE(rtmp); 163 ierr = ISRestoreIndices(isicol,&ic); CHKERR(ierr); 164 ierr = ISRestoreIndices(isrow,&r); CHKERR(ierr); 165 ierr = ISDestroy(isicol); CHKERR(ierr); 166 fact->factor = FACTOR_LU; 167 168 return 0; 169 } 170 int MatiAIJSolve(Mat mat,Vec bb, Vec xx) 171 { 172 Matiaij *aij = (Matiaij *) mat->data; 173 IS iscol = mat->col, isrow = mat->row; 174 int *r,*c, ierr, i, j, n = aij->m, *vi, *ai = aij->i, *aj = aij->j; 175 int nz; 176 Scalar *x,*b,*tmp, *aa = aij->a, sum, *v; 177 178 if (ierr = VecGetArray(bb,&b)) SETERR(ierr,0); 179 if (ierr = VecGetArray(xx,&x)) SETERR(ierr,0); 180 tmp = (Scalar *) MALLOC(n*sizeof(Scalar)); CHKPTR(tmp); 181 182 if (ierr = ISGetIndices(isrow,&r)) SETERR(ierr,0); 183 if (ierr = ISGetIndices(iscol,&c)) SETERR(ierr,0); c = c + (n-1); 184 185 /* forward solve the lower triangular */ 186 tmp[0] = b[*r++]; 187 for ( i=1; i<n; i++ ) { 188 v = aa + ai[i] - 1; 189 vi = aj + ai[i] - 1; 190 nz = aij->diag[i] - ai[i]; 191 sum = b[*r++]; 192 while (nz--) sum -= *v++ * tmp[*vi++ - 1]; 193 tmp[i] = sum; 194 } 195 196 /* backward solve the upper triangular */ 197 for ( i=n-1; i>=0; i-- ){ 198 v = aa + aij->diag[i]; 199 vi = aj + aij->diag[i]; 200 nz = ai[i+1] - aij->diag[i] - 1; 201 sum = tmp[i]; 202 while (nz--) sum -= *v++ * tmp[*vi++ - 1]; 203 x[*c--] = tmp[i] = sum*aa[aij->diag[i]-1]; 204 } 205 206 FREE(tmp); 207 return 0; 208 } 209