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