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