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