1 2 3 #include <petsc.h> 4 #include <petscksp.h> 5 #include "private/kspimpl.h" 6 #include "petscpc.h" 7 #include "../src/mat/impls/aij/seq/aij.h" 8 #include "../src/mat/impls/sbaij/seq/sbaij.h" 9 #include "../src/mat/impls/aij/seq/bas/spbas.h" 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatICCFactorSymbolic_SeqAIJ_Bas" 13 PetscErrorCode MatICCFactorSymbolic_SeqAIJ_Bas(Mat fact,Mat A,IS perm,const MatFactorInfo *info) 14 { 15 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 16 Mat_SeqSBAIJ *b; 17 PetscErrorCode ierr; 18 PetscTruth perm_identity,missing; 19 PetscInt reallocs=0,i,*ai=a->i,*aj=a->j,am=A->rmap->n,*ui; 20 const PetscInt *rip,*riip; 21 PetscInt j; 22 PetscInt d; 23 PetscInt ncols,*cols,*uj; 24 PetscReal fill=info->fill,levels=info->levels; 25 IS iperm; 26 spbas_matrix Pattern_0, Pattern_P; 27 28 PetscFunctionBegin; 29 if (A->rmap->n != A->cmap->n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",A->rmap->n,A->cmap->n); 30 ierr = MatMissingDiagonal(A,&missing,&d);CHKERRQ(ierr); 31 if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",d); 32 ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr); 33 ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr); 34 35 36 /* ICC(0) without matrix ordering: simply copies fill pattern */ 37 if (!levels && perm_identity) { 38 ierr = PetscMalloc((am+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr); 39 ui[0] = 0; 40 41 for (i=0; i<am; i++) { 42 ui[i+1] = ui[i] + ai[i+1] - a->diag[i]; 43 } 44 ierr = PetscMalloc((ui[am]+1)*sizeof(PetscInt),&uj);CHKERRQ(ierr); 45 cols = uj; 46 for (i=0; i<am; i++) { 47 aj = a->j + a->diag[i]; 48 ncols = ui[i+1] - ui[i]; 49 for (j=0; j<ncols; j++) *cols++ = *aj++; 50 } 51 } else { /* case: levels>0 || (levels=0 && !perm_identity) */ 52 ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr); 53 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr); 54 55 // Create spbas_matrix for pattern 56 ierr = spbas_pattern_only(am, am, ai, aj, &Pattern_0); CHKERRQ(ierr); 57 58 // Apply the permutation 59 ierr = spbas_apply_reordering( &Pattern_0, rip, riip); CHKERRQ(ierr); 60 61 // Raise the power 62 ierr = spbas_power( Pattern_0, (int) levels+1, &Pattern_P); 63 CHKERRQ(ierr); 64 ierr = spbas_delete( Pattern_0 ); CHKERRQ(ierr); 65 66 // Keep only upper triangle of pattern 67 ierr = spbas_keep_upper( &Pattern_P ); 68 69 // Convert to Sparse Row Storage 70 ierr = spbas_matrix_to_crs(Pattern_P, NULL, &ui, &uj); CHKERRQ(ierr); 71 ierr = spbas_delete(Pattern_P);CHKERRQ(ierr); 72 } /* end of case: levels>0 || (levels=0 && !perm_identity) */ 73 74 /* put together the new matrix in MATSEQSBAIJ format */ 75 76 b = (Mat_SeqSBAIJ*)(fact)->data; 77 b->singlemalloc = PETSC_FALSE; 78 ierr = PetscMalloc((ui[am]+1)*sizeof(MatScalar),&b->a);CHKERRQ(ierr); 79 b->j = uj; 80 b->i = ui; 81 b->diag = 0; 82 b->ilen = 0; 83 b->imax = 0; 84 b->row = perm; 85 b->col = perm; 86 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 87 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 88 b->icol = iperm; 89 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 90 ierr = PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr); 91 ierr = PetscLogObjectMemory((fact),(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 92 b->maxnz = b->nz = ui[am]; 93 b->free_a = PETSC_TRUE; 94 b->free_ij = PETSC_TRUE; 95 96 (fact)->info.factor_mallocs = reallocs; 97 (fact)->info.fill_ratio_given = fill; 98 if (ai[am] != 0) { 99 (fact)->info.fill_ratio_needed = ((PetscReal)ui[am])/((PetscReal)ai[am]); 100 } else { 101 (fact)->info.fill_ratio_needed = 0.0; 102 } 103 /* (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_inplace; */ 104 PetscFunctionReturn(0); 105 } 106 107 108 #undef __FUNCT__ 109 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqAIJ_Bas" 110 PetscErrorCode MatCholeskyFactorNumeric_SeqAIJ_Bas(Mat B,Mat A,const MatFactorInfo *info) 111 { 112 Mat C = B; 113 Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data; 114 IS ip=b->row,iip = b->icol; 115 PetscErrorCode ierr; 116 const PetscInt *rip,*riip; 117 PetscInt mbs=A->rmap->n,*bi=b->i,*bj=b->j; 118 119 MatScalar *ba=b->a; 120 PetscReal shiftnz = info->shiftnz; 121 PetscReal droptol = -1; 122 PetscTruth perm_identity; 123 spbas_matrix Pattern, matrix_L,matrix_LT; 124 PetscReal mem_reduction; 125 126 PetscFunctionBegin; 127 // Reduce memory requirements: 128 // erase values of B-matrix 129 ierr = PetscFree(ba); CHKERRQ(ierr); 130 // Compress (maximum) sparseness pattern of B-matrix 131 ierr = spbas_compress_pattern(bi, bj, mbs, mbs, SPBAS_DIAGONAL_OFFSETS,&Pattern, &mem_reduction);CHKERRQ(ierr); 132 ierr = PetscFree(bi); CHKERRQ(ierr); 133 ierr = PetscFree(bj); CHKERRQ(ierr); 134 135 ierr = PetscInfo1(PETSC_NULL," compression rate for spbas_compress_pattern %G \n",mem_reduction);CHKERRQ(ierr); 136 137 // Make Cholesky decompositions with larger Manteuffel shifts until no more 138 // negative diagonals are found. 139 ierr = ISGetIndices(ip,&rip);CHKERRQ(ierr); 140 ierr = ISGetIndices(iip,&riip);CHKERRQ(ierr); 141 142 if (info->usedt) { 143 droptol = info->dt; 144 } 145 for (ierr = NEGATIVE_DIAGONAL; ierr == NEGATIVE_DIAGONAL; ) 146 { 147 ierr = spbas_incomplete_cholesky( A, rip, riip, Pattern, droptol, shiftnz, 148 &matrix_LT); 149 if (ierr == NEGATIVE_DIAGONAL) 150 { 151 shiftnz *= 1.5; 152 if (shiftnz < 1e-5) shiftnz=1e-5; 153 ierr = PetscInfo1(PETSC_NULL,"spbas_incomplete_cholesky found a negative diagonal. Trying again with Manteuffel shift=%G\n",shiftnz);CHKERRQ(ierr); 154 } 155 } 156 CHKERRQ(ierr); 157 ierr = spbas_delete(Pattern); CHKERRQ(ierr); 158 159 ierr = PetscInfo1(PETSC_NULL," memory_usage for spbas_incomplete_cholesky %G bytes per row\n", 160 (PetscReal) spbas_memory_requirement( matrix_LT)/ (PetscReal) mbs);CHKERRQ(ierr); 161 162 ierr = ISRestoreIndices(ip,&rip);CHKERRQ(ierr); 163 ierr = ISRestoreIndices(iip,&riip);CHKERRQ(ierr); 164 165 // Convert spbas_matrix to compressed row storage 166 ierr = spbas_transpose(matrix_LT, &matrix_L); CHKERRQ(ierr); 167 ierr = spbas_delete(matrix_LT); CHKERRQ(ierr); 168 #if defined(foo) 169 { ierr = spbas_dump("factorL",matrix_L); CHKERRQ(ierr);} 170 #endif 171 ierr = spbas_matrix_to_crs(matrix_L, &ba, &bi, &bj); CHKERRQ(ierr); 172 b->i=bi; b->j=bj; b->a=ba; 173 ierr = spbas_delete(matrix_L); CHKERRQ(ierr); 174 175 // Set the appropriate solution functions 176 ierr = ISIdentity(ip,&perm_identity);CHKERRQ(ierr); 177 if (perm_identity){ 178 (B)->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 179 (B)->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 180 (B)->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 181 (B)->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 182 } else { 183 (B)->ops->solve = MatSolve_SeqSBAIJ_1_inplace; 184 (B)->ops->solvetranspose = MatSolve_SeqSBAIJ_1_inplace; 185 (B)->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_1_inplace; 186 (B)->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_1_inplace; 187 } 188 189 C->assembled = PETSC_TRUE; 190 C->preallocated = PETSC_TRUE; 191 ierr = PetscLogFlops(C->rmap->n);CHKERRQ(ierr); 192 193 // Optionally, print the factor matrix to file 194 #if defined(foo) 195 { 196 FILE * factfile = fopen("factorL","w"); 197 if (!factfile) CHKERRQ((PetscErrorCode) 10); 198 for (i=0; i<mbs; i++) 199 { 200 for (j=bi[i]; j<bi[i+1]; j++) 201 { 202 fprintf(factfile,"%d %d %e\n",i,bj[j],ba[j]); 203 } 204 } 205 fclose(factfile); 206 } 207 #endif 208 PetscFunctionReturn(0); 209 } 210 211 EXTERN_C_BEGIN 212 #undef __FUNCT__ 213 #define __FUNCT__ "MatGetFactor_seqaij_bas" 214 PetscErrorCode MatGetFactor_seqaij_bas(Mat A,MatFactorType ftype,Mat *B) 215 { 216 PetscInt n = A->rmap->n; 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 221 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 222 if (ftype == MAT_FACTOR_ICC) { 223 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 224 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 225 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqAIJ_Bas; 226 (*B)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJ_Bas; 227 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 228 (*B)->factor = ftype; 229 PetscFunctionReturn(0); 230 } 231 EXTERN_C_END 232