1 2 3 #include "aij.h" 4 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; 101 aijnew->a = (Scalar *) MALLOC( len ); CHKPTR(aijnew->a); 102 aijnew->j = ajnew; 103 aijnew->i = ainew; 104 aij->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,m,fm, nnz, nzi,len; 118 Scalar *rtmp,*vnew,*v; 119 120 if (ierr = ISInvertPermutation(iscol,&isicol)) SETERR(ierr,0); 121 ISGetIndices(isrow,&r); ISGetIndices(isicol,&ic); 122 rtmp = (Scalar *) MALLOC( (n+1)*sizeof(Scalar) ); CHKPTR(rtmp); 123 124 for ( i=0; i<n; i++ ) { 125 nz = ai[i+1] - ai[i]; 126 ajtmp = aj + ai[i] - 1; 127 vnew = aij->a + 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[i+1] - aij->i[i]; 132 ajtmpold = aij->j + aij->i[i] - 1; 133 v = aij->a + aij->i[i] - 1; 134 for ( j=0; j<nz; j++ ) rtmp[ic[ajtmpold[j]-1]] = v[i]; 135 136 row = *ajtmp++: 137 while (row < i) { 138 pc = rm 139 } 140 } 141 ISDestroy(isicol); 142 143 return 0; 144 } 145