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