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 int *r,*ic,ierr,i,n = a->mbs,*ai = a->i,*aj = a->j; 19 int *ainew,*ajnew,jmax,*fill,*ajtmp,nz,bs = a->bs,bs2=a->bs2; 20 int *idnew,idx,row,m,fm,nnz,nzi,realloc = 0,nzbd,*im; 21 PetscReal f = 1.0; 22 23 PetscFunctionBegin; 24 if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 25 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 26 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 27 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 28 29 f = info->fill; 30 /* get new row pointers */ 31 ierr = PetscMalloc((n+1)*sizeof(int),&ainew);CHKERRQ(ierr); 32 ainew[0] = 0; 33 /* don't know how many column pointers are needed so estimate */ 34 jmax = (int)(f*ai[n] + 1); 35 ierr = PetscMalloc((jmax)*sizeof(int),&ajnew);CHKERRQ(ierr); 36 /* fill is a linked list of nonzeros in active row */ 37 ierr = PetscMalloc((2*n+1)*sizeof(int),&fill);CHKERRQ(ierr); 38 im = fill + n + 1; 39 /* idnew is location of diagonal in factor */ 40 ierr = PetscMalloc((n+1)*sizeof(int),&idnew);CHKERRQ(ierr); 41 idnew[0] = 0; 42 43 for (i=0; i<n; i++) { 44 /* first copy previous fill into linked list */ 45 nnz = nz = ai[r[i]+1] - ai[r[i]]; 46 if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix"); 47 ajtmp = aj + ai[r[i]]; 48 fill[n] = n; 49 while (nz--) { 50 fm = n; 51 idx = ic[*ajtmp++]; 52 do { 53 m = fm; 54 fm = fill[m]; 55 } while (fm < idx); 56 fill[m] = idx; 57 fill[idx] = fm; 58 } 59 row = fill[n]; 60 while (row < i) { 61 ajtmp = ajnew + idnew[row] + 1; 62 nzbd = 1 + idnew[row] - ainew[row]; 63 nz = im[row] - nzbd; 64 fm = row; 65 while (nz-- > 0) { 66 idx = *ajtmp++; 67 nzbd++; 68 if (idx == i) im[row] = nzbd; 69 do { 70 m = fm; 71 fm = fill[m]; 72 } while (fm < idx); 73 if (fm != idx) { 74 fill[m] = idx; 75 fill[idx] = fm; 76 fm = idx; 77 nnz++; 78 } 79 } 80 row = fill[row]; 81 } 82 /* copy new filled row into permanent storage */ 83 ainew[i+1] = ainew[i] + nnz; 84 if (ainew[i+1] > jmax) { 85 86 /* estimate how much additional space we will need */ 87 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 88 /* just double the memory each time */ 89 int maxadd = jmax; 90 /* maxadd = (int)((f*(ai[n]+1)*(n-i+5))/n); */ 91 if (maxadd < nnz) maxadd = (n-i)*(nnz+1); 92 jmax += maxadd; 93 94 /* allocate a longer ajnew */ 95 ierr = PetscMalloc(jmax*sizeof(int),&ajtmp);CHKERRQ(ierr); 96 ierr = PetscMemcpy(ajtmp,ajnew,ainew[i]*sizeof(int));CHKERRQ(ierr); 97 ierr = PetscFree(ajnew);CHKERRQ(ierr); 98 ajnew = ajtmp; 99 realloc++; /* count how many times we realloc */ 100 } 101 ajtmp = ajnew + ainew[i]; 102 fm = fill[n]; 103 nzi = 0; 104 im[i] = nnz; 105 while (nnz--) { 106 if (fm < i) nzi++; 107 *ajtmp++ = fm; 108 fm = fill[fm]; 109 } 110 idnew[i] = ainew[i] + nzi; 111 } 112 113 if (ai[n] != 0) { 114 PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]); 115 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Reallocs %d Fill ratio:given %g needed %g\n",realloc,f,af); 116 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Run with -pc_lu_fill %g or use \n",af); 117 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:PCLUSetFill(pc,%g);\n",af); 118 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:for best performance.\n"); 119 } else { 120 PetscLogInfo(A,"MatLUFactorSymbolic_SeqBAIJ:Empty matrix.\n"); 121 } 122 123 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 124 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 125 126 ierr = PetscFree(fill);CHKERRQ(ierr); 127 128 /* put together the new matrix */ 129 ierr = MatCreate(A->comm,bs*n,bs*n,bs*n,bs*n,B);CHKERRQ(ierr); 130 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 131 ierr = MatSeqBAIJSetPreallocation(*B,bs,0,PETSC_NULL);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