1 /*$Id: baijfact3.c,v 1.6 2001/08/06 21:15:36 bsmith Exp $*/ 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 __FUNCT__ 15 #define __FUNCT__ "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 if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 27 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 28 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 29 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 30 31 if (info) f = info->fill; 32 /* get new row pointers */ 33 ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr); 34 ainew[0] = 0; 35 /* don't know how many column pointers are needed so estimate */ 36 jmax = (int)(f*ai[n] + 1); 37 ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr); 38 /* fill is a linked list of nonzeros in active row */ 39 ierr = PetscMalloc((2*n+1)*sizeof(int),&fill);CHKERRQ(ierr); 40 im = fill + n + 1; 41 /* idnew is location of diagonal in factor */ 42 ierr = PetscMalloc((n+1)*sizeof(int),&idnew);CHKERRQ(ierr); 43 idnew[0] = 0; 44 45 for (i=0; i<n; i++) { 46 /* first copy previous fill into linked list */ 47 nnz = nz = ai[r[i]+1] - ai[r[i]]; 48 if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix"); 49 ajtmp = aj + ai[r[i]]; 50 fill[n] = n; 51 while (nz--) { 52 fm = n; 53 idx = ic[*ajtmp++]; 54 do { 55 m = fm; 56 fm = fill[m]; 57 } while (fm < idx); 58 fill[m] = idx; 59 fill[idx] = fm; 60 } 61 row = fill[n]; 62 while (row < i) { 63 ajtmp = ajnew + idnew[row] + 1; 64 nzbd = 1 + idnew[row] - ainew[row]; 65 nz = im[row] - nzbd; 66 fm = row; 67 while (nz-- > 0) { 68 idx = *ajtmp++; 69 nzbd++; 70 if (idx == i) im[row] = nzbd; 71 do { 72 m = fm; 73 fm = fill[m]; 74 } while (fm < idx); 75 if (fm != idx) { 76 fill[m] = idx; 77 fill[idx] = fm; 78 fm = idx; 79 nnz++; 80 } 81 } 82 row = fill[row]; 83 } 84 /* copy new filled row into permanent storage */ 85 ainew[i+1] = ainew[i] + nnz; 86 if (ainew[i+1] > jmax) { 87 88 /* estimate how much additional space we will need */ 89 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 90 /* just double the memory each time */ 91 int maxadd = jmax; 92 /* maxadd = (int)((f*(ai[n]+1)*(n-i+5))/n); */ 93 if (maxadd < nnz) maxadd = (n-i)*(nnz+1); 94 jmax += maxadd; 95 96 /* allocate a longer ajnew */ 97 ierr = PetscMalloc(jmax*sizeof(int),&ajtmp);CHKERRQ(ierr); 98 ierr = PetscMemcpy(ajtmp,ajnew,ainew[i]*sizeof(int));CHKERRQ(ierr); 99 ierr = PetscFree(ajnew);CHKERRQ(ierr); 100 ajnew = ajtmp; 101 realloc++; /* count how many times we realloc */ 102 } 103 ajtmp = ajnew + ainew[i]; 104 fm = fill[n]; 105 nzi = 0; 106 im[i] = nnz; 107 while (nnz--) { 108 if (fm < i) nzi++; 109 *ajtmp++ = fm; 110 fm = fill[fm]; 111 } 112 idnew[i] = ainew[i] + nzi; 113 } 114 115 if (ai[n] != 0) { 116 PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]); 117 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af); 118 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Run with -pc_lu_fill %g or use \n",af); 119 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:PCLUSetFill(pc,%g);\n",af); 120 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:for best performance.\n"); 121 } else { 122 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Empty matrix.\n"); 123 } 124 125 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 126 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 127 128 ierr = PetscFree(fill);CHKERRQ(ierr); 129 130 /* put together the new matrix */ 131 ierr = MatCreateSeqBAIJ(A->comm,bs,bs*n,bs*n,0,PETSC_NULL,B);CHKERRQ(ierr); 132 PetscLogObjectParent(*B,isicol); 133 b = (Mat_SeqBAIJ*)(*B)->data; 134 ierr = PetscFree(b->imax);CHKERRQ(ierr); 135 b->singlemalloc = PETSC_FALSE; 136 /* the next line frees the default space generated by the Create() */ 137 ierr = PetscFree(b->a);CHKERRQ(ierr); 138 ierr = PetscFree(b->ilen);CHKERRQ(ierr); 139 ierr = PetscMalloc((ainew[n]+1)*sizeof(MatScalar)*bs2,&b->a);CHKERRQ(ierr); 140 b->j = ajnew; 141 b->i = ainew; 142 b->diag = idnew; 143 b->ilen = 0; 144 b->imax = 0; 145 b->row = isrow; 146 b->col = iscol; 147 b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE; 148 ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr); 149 ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr); 150 b->icol = isicol; 151 ierr = PetscMalloc((bs*n+bs)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr); 152 /* In b structure: Free imax, ilen, old a, old j. 153 Allocate idnew, solve_work, new a, new j */ 154 PetscLogObjectMemory(*B,(ainew[n]-n)*(sizeof(int)+sizeof(MatScalar))); 155 b->maxnz = b->nz = ainew[n]; 156 157 (*B)->factor = FACTOR_LU; 158 (*B)->info.factor_mallocs = realloc; 159 (*B)->info.fill_ratio_given = f; 160 if (ai[n] != 0) { 161 (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]); 162 } else { 163 (*B)->info.fill_ratio_needed = 0.0; 164 } 165 PetscFunctionReturn(0); 166 } 167