1 /*$Id: baijfact3.c,v 1.2 2001/01/15 21:45:50 bsmith Exp balay $*/ 2 /* 3 Factorization code for BAIJ format. 4 */ 5 #include "src/mat/impls/baij/seq/baij.h" 6 #include "src/vec/vecimpl.h" 7 #include "src/inline/ilu.h" 8 9 /* 10 The symbolic factorization code is identical to that for AIJ format, 11 except for very small changes since this is now a SeqBAIJ datastructure. 12 NOT good code reuse. 13 */ 14 #undef __FUNC__ 15 #define __FUNC__ "MatLUFactorSymbolic_SeqBAIJ" 16 int MatLUFactorSymbolic_SeqBAIJ(Mat A,IS isrow,IS iscol,MatLUInfo *info,Mat *B) 17 { 18 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b; 19 IS isicol; 20 int *r,*ic,ierr,i,n = a->mbs,*ai = a->i,*aj = a->j; 21 int *ainew,*ajnew,jmax,*fill,*ajtmp,nz,bs = a->bs,bs2=a->bs2; 22 int *idnew,idx,row,m,fm,nnz,nzi,realloc = 0,nzbd,*im; 23 PetscReal f = 1.0; 24 25 PetscFunctionBegin; 26 PetscValidHeaderSpecific(isrow,IS_COOKIE); 27 PetscValidHeaderSpecific(iscol,IS_COOKIE); 28 if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 29 30 if (!isrow) { 31 ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&isrow);CHKERRQ(ierr); 32 } 33 if (!iscol) { 34 ierr = ISCreateStride(PETSC_COMM_SELF,A->M,0,1,&iscol);CHKERRQ(ierr); 35 } 36 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 37 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 38 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 39 40 if (info) f = info->fill; 41 /* get new row pointers */ 42 ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr); 43 ainew[0] = 0; 44 /* don't know how many column pointers are needed so estimate */ 45 jmax = (int)(f*ai[n] + 1); 46 ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr); 47 /* fill is a linked list of nonzeros in active row */ 48 ierr = PetscMalloc((2*n+1)*sizeof(int),&fill);CHKERRQ(ierr); 49 im = fill + n + 1; 50 /* idnew is location of diagonal in factor */ 51 ierr = PetscMalloc((n+1)*sizeof(int),&idnew);CHKERRQ(ierr); 52 idnew[0] = 0; 53 54 for (i=0; i<n; i++) { 55 /* first copy previous fill into linked list */ 56 nnz = nz = ai[r[i]+1] - ai[r[i]]; 57 if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix"); 58 ajtmp = aj + ai[r[i]]; 59 fill[n] = n; 60 while (nz--) { 61 fm = n; 62 idx = ic[*ajtmp++]; 63 do { 64 m = fm; 65 fm = fill[m]; 66 } while (fm < idx); 67 fill[m] = idx; 68 fill[idx] = fm; 69 } 70 row = fill[n]; 71 while (row < i) { 72 ajtmp = ajnew + idnew[row] + 1; 73 nzbd = 1 + idnew[row] - ainew[row]; 74 nz = im[row] - nzbd; 75 fm = row; 76 while (nz-- > 0) { 77 idx = *ajtmp++; 78 nzbd++; 79 if (idx == i) im[row] = nzbd; 80 do { 81 m = fm; 82 fm = fill[m]; 83 } while (fm < idx); 84 if (fm != idx) { 85 fill[m] = idx; 86 fill[idx] = fm; 87 fm = idx; 88 nnz++; 89 } 90 } 91 row = fill[row]; 92 } 93 /* copy new filled row into permanent storage */ 94 ainew[i+1] = ainew[i] + nnz; 95 if (ainew[i+1] > jmax) { 96 97 /* estimate how much additional space we will need */ 98 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 99 /* just double the memory each time */ 100 int maxadd = jmax; 101 /* maxadd = (int)((f*(ai[n]+1)*(n-i+5))/n); */ 102 if (maxadd < nnz) maxadd = (n-i)*(nnz+1); 103 jmax += maxadd; 104 105 /* allocate a longer ajnew */ 106 ierr = PetscMalloc(jmax*sizeof(int),&ajtmp);CHKERRQ(ierr); 107 ierr = PetscMemcpy(ajtmp,ajnew,ainew[i]*sizeof(int));CHKERRQ(ierr); 108 ierr = PetscFree(ajnew);CHKERRQ(ierr); 109 ajnew = ajtmp; 110 realloc++; /* count how many times we realloc */ 111 } 112 ajtmp = ajnew + ainew[i]; 113 fm = fill[n]; 114 nzi = 0; 115 im[i] = nnz; 116 while (nnz--) { 117 if (fm < i) nzi++; 118 *ajtmp++ = fm; 119 fm = fill[fm]; 120 } 121 idnew[i] = ainew[i] + nzi; 122 } 123 124 if (ai[n] != 0) { 125 PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]); 126 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af); 127 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Run with -pc_lu_fill %g or use \n",af); 128 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:PCLUSetFill(pc,%g);\n",af); 129 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:for best performance.\n"); 130 } else { 131 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Empty matrix.\n"); 132 } 133 134 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 135 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 136 137 ierr = PetscFree(fill);CHKERRQ(ierr); 138 139 /* put together the new matrix */ 140 ierr = MatCreateSeqBAIJ(A->comm,bs,bs*n,bs*n,0,PETSC_NULL,B);CHKERRQ(ierr); 141 PetscLogObjectParent(*B,isicol); 142 b = (Mat_SeqBAIJ*)(*B)->data; 143 ierr = PetscFree(b->imax);CHKERRQ(ierr); 144 b->singlemalloc = PETSC_FALSE; 145 /* the next line frees the default space generated by the Create() */ 146 ierr = PetscFree(b->a);CHKERRQ(ierr); 147 ierr = PetscFree(b->ilen);CHKERRQ(ierr); 148 ierr = PetscMalloc((ainew[n]+1)*sizeof(MatScalar)*bs2,&b->a);CHKERRQ(ierr); 149 b->j = ajnew; 150 b->i = ainew; 151 b->diag = idnew; 152 b->ilen = 0; 153 b->imax = 0; 154 b->row = isrow; 155 b->col = iscol; 156 ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr); 157 ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr); 158 b->icol = isicol; 159 ierr = PetscMalloc((bs*n+bs)*sizeof(Scalar),&b->solve_work);CHKERRQ(ierr); 160 /* In b structure: Free imax, ilen, old a, old j. 161 Allocate idnew, solve_work, new a, new j */ 162 PetscLogObjectMemory(*B,(ainew[n]-n)*(sizeof(int)+sizeof(MatScalar))); 163 b->maxnz = b->nz = ainew[n]; 164 165 (*B)->factor = FACTOR_LU; 166 (*B)->info.factor_mallocs = realloc; 167 (*B)->info.fill_ratio_given = f; 168 if (ai[n] != 0) { 169 (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]); 170 } else { 171 (*B)->info.fill_ratio_needed = 0.0; 172 } 173 PetscFunctionReturn(0); 174 } 175