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