1 2 /* 3 Provides an interface to the CHOLMOD sparse solver available through SuiteSparse version 4.2.1 4 5 When built with PETSC_USE_64BIT_INDICES this will use Suitesparse_long as the 6 integer type in UMFPACK, otherwise it will use int. This means 7 all integers in this file as simply declared as PetscInt. Also it means 8 that one cannot use 64BIT_INDICES on 32bit machines [as Suitesparse_long is 32bit only] 9 10 */ 11 12 #include <../src/mat/impls/sbaij/seq/sbaij.h> 13 #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h> 14 15 /* 16 This is a terrible hack, but it allows the error handler to retain a context. 17 Note that this hack really cannot be made both reentrant and concurrent. 18 */ 19 static Mat static_F; 20 21 static void CholmodErrorHandler(int status,const char *file,int line,const char *message) 22 { 23 PetscFunctionBegin; 24 if (status > CHOLMOD_OK) { 25 PetscCallVoid(PetscInfo(static_F,"CHOLMOD warning %d at %s:%d: %s\n",status,file,line,message)); 26 } else if (status == CHOLMOD_OK) { /* Documentation says this can happen, but why? */ 27 PetscCallVoid(PetscInfo(static_F,"CHOLMOD OK at %s:%d: %s\n",file,line,message)); 28 } else { 29 PetscCallVoid(PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n",status,file,line,message)); 30 } 31 PetscFunctionReturnVoid(); 32 } 33 34 PetscErrorCode CholmodStart(Mat F) 35 { 36 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 37 cholmod_common *c; 38 PetscBool flg; 39 40 PetscFunctionBegin; 41 if (chol->common) PetscFunctionReturn(0); 42 PetscCall(PetscMalloc1(1,&chol->common)); 43 PetscCall(!cholmod_X_start(chol->common)); 44 45 c = chol->common; 46 c->error_handler = CholmodErrorHandler; 47 48 #define CHOLMOD_OPTION_DOUBLE(name,help) do { \ 49 PetscReal tmp = (PetscReal)c->name; \ 50 PetscCall(PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 51 c->name = (double)tmp; \ 52 } while (0) 53 54 #define CHOLMOD_OPTION_INT(name,help) do { \ 55 PetscInt tmp = (PetscInt)c->name; \ 56 PetscCall(PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 57 c->name = (int)tmp; \ 58 } while (0) 59 60 #define CHOLMOD_OPTION_SIZE_T(name,help) do { \ 61 PetscReal tmp = (PetscInt)c->name; \ 62 PetscCall(PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 63 PetscCheck(tmp >= 0,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"value must be positive"); \ 64 c->name = (size_t)tmp; \ 65 } while (0) 66 67 #define CHOLMOD_OPTION_BOOL(name,help) do { \ 68 PetscBool tmp = (PetscBool) !!c->name; \ 69 PetscCall(PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 70 c->name = (int)tmp; \ 71 } while (0) 72 73 PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat"); 74 CHOLMOD_OPTION_INT(nmethods,"Number of different ordering methods to try"); 75 76 #if defined(PETSC_USE_SUITESPARSE_GPU) 77 c->useGPU = 1; 78 CHOLMOD_OPTION_INT(useGPU,"Use GPU for BLAS 1, otherwise 0"); 79 CHOLMOD_OPTION_SIZE_T(maxGpuMemBytes,"Maximum memory to allocate on the GPU"); 80 CHOLMOD_OPTION_DOUBLE(maxGpuMemFraction,"Fraction of available GPU memory to allocate"); 81 #endif 82 83 /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */ 84 chol->pack = (PetscBool)c->final_pack; 85 PetscCall(PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL)); 86 c->final_pack = (int)chol->pack; 87 88 CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D"); 89 CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified"); 90 CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified"); 91 CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified"); 92 CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]"); 93 { 94 static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0}; 95 PetscCall(PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL)); 96 } 97 if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization"); 98 CHOLMOD_OPTION_BOOL(final_asis,"Leave factors \"as is\""); 99 CHOLMOD_OPTION_BOOL(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)"); 100 if (!c->final_asis) { 101 CHOLMOD_OPTION_BOOL(final_super,"Leave supernodal factors instead of converting to simplicial"); 102 CHOLMOD_OPTION_BOOL(final_ll,"Turn LDL' factorization into LL'"); 103 CHOLMOD_OPTION_BOOL(final_monotonic,"Ensure columns are monotonic when done"); 104 CHOLMOD_OPTION_BOOL(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation"); 105 } 106 { 107 PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]}; 108 PetscInt n = 3; 109 PetscCall(PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg)); 110 PetscCheck(!flg || n == 3,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_zrelax"); 111 if (flg) while (n--) c->zrelax[n] = (double)tmp[n]; 112 } 113 { 114 PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]}; 115 PetscCall(PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg)); 116 PetscCheck(!flg || n == 3,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_nrelax"); 117 if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n]; 118 } 119 CHOLMOD_OPTION_BOOL(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]"); 120 CHOLMOD_OPTION_BOOL(default_nesdis,"Use NESDIS instead of METIS for nested dissection"); 121 CHOLMOD_OPTION_INT(print,"Verbosity level"); 122 PetscOptionsEnd(); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc,PetscBool *valloc) 127 { 128 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)A->data; 129 PetscBool vallocin = PETSC_FALSE; 130 131 PetscFunctionBegin; 132 PetscCall(PetscMemzero(C,sizeof(*C))); 133 /* CHOLMOD uses column alignment, SBAIJ stores the upper factor, so we pass it on as a lower factor, swapping the meaning of row and column */ 134 C->nrow = (size_t)A->cmap->n; 135 C->ncol = (size_t)A->rmap->n; 136 C->nzmax = (size_t)sbaij->maxnz; 137 C->p = sbaij->i; 138 C->i = sbaij->j; 139 if (values) { 140 #if defined(PETSC_USE_COMPLEX) 141 /* we need to pass CHOLMOD the conjugate matrix */ 142 PetscScalar *v; 143 PetscInt i; 144 145 PetscCall(PetscMalloc1(sbaij->maxnz,&v)); 146 for (i = 0; i < sbaij->maxnz; i++) v[i] = PetscConj(sbaij->a[i]); 147 C->x = v; 148 vallocin = PETSC_TRUE; 149 #else 150 C->x = sbaij->a; 151 #endif 152 } 153 C->stype = -1; 154 C->itype = CHOLMOD_INT_TYPE; 155 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 156 C->dtype = CHOLMOD_DOUBLE; 157 C->sorted = 1; 158 C->packed = 1; 159 *aijalloc = PETSC_FALSE; 160 *valloc = vallocin; 161 PetscFunctionReturn(0); 162 } 163 164 #define GET_ARRAY_READ 0 165 #define GET_ARRAY_WRITE 1 166 167 PetscErrorCode VecWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y) 168 { 169 PetscScalar *x; 170 PetscInt n; 171 172 PetscFunctionBegin; 173 PetscCall(PetscMemzero(Y,sizeof(*Y))); 174 switch (rw) { 175 case GET_ARRAY_READ: 176 PetscCall(VecGetArrayRead(X,(const PetscScalar**)&x)); 177 break; 178 case GET_ARRAY_WRITE: 179 PetscCall(VecGetArrayWrite(X,&x)); 180 break; 181 default: 182 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 183 break; 184 } 185 PetscCall(VecGetSize(X,&n)); 186 187 Y->x = x; 188 Y->nrow = n; 189 Y->ncol = 1; 190 Y->nzmax = n; 191 Y->d = n; 192 Y->xtype = CHOLMOD_SCALAR_TYPE; 193 Y->dtype = CHOLMOD_DOUBLE; 194 PetscFunctionReturn(0); 195 } 196 197 PetscErrorCode VecUnWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y) 198 { 199 PetscFunctionBegin; 200 switch (rw) { 201 case GET_ARRAY_READ: 202 PetscCall(VecRestoreArrayRead(X,(const PetscScalar**)&Y->x)); 203 break; 204 case GET_ARRAY_WRITE: 205 PetscCall(VecRestoreArrayWrite(X,(PetscScalar**)&Y->x)); 206 break; 207 default: 208 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 209 break; 210 } 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode MatDenseWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 215 { 216 PetscScalar *x; 217 PetscInt m,n,lda; 218 219 PetscFunctionBegin; 220 PetscCall(PetscMemzero(Y,sizeof(*Y))); 221 switch (rw) { 222 case GET_ARRAY_READ: 223 PetscCall(MatDenseGetArrayRead(X,(const PetscScalar**)&x)); 224 break; 225 case GET_ARRAY_WRITE: 226 PetscCall(MatDenseGetArrayWrite(X,&x)); 227 break; 228 default: 229 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 230 break; 231 } 232 PetscCall(MatDenseGetLDA(X,&lda)); 233 PetscCall(MatGetLocalSize(X,&m,&n)); 234 235 Y->x = x; 236 Y->nrow = m; 237 Y->ncol = n; 238 Y->nzmax = lda*n; 239 Y->d = lda; 240 Y->xtype = CHOLMOD_SCALAR_TYPE; 241 Y->dtype = CHOLMOD_DOUBLE; 242 PetscFunctionReturn(0); 243 } 244 245 PetscErrorCode MatDenseUnWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 246 { 247 PetscFunctionBegin; 248 switch (rw) { 249 case GET_ARRAY_READ: 250 PetscCall(MatDenseRestoreArrayRead(X,(const PetscScalar**)&Y->x)); 251 break; 252 case GET_ARRAY_WRITE: 253 /* we don't have MatDenseRestoreArrayWrite */ 254 PetscCall(MatDenseRestoreArray(X,(PetscScalar**)&Y->x)); 255 break; 256 default: 257 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 258 break; 259 } 260 PetscFunctionReturn(0); 261 } 262 263 PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F) 264 { 265 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 266 267 PetscFunctionBegin; 268 if (chol->spqrfact) { 269 PetscCall(!SuiteSparseQR_C_free(&chol->spqrfact, chol->common)); 270 } 271 if (chol->factor) { 272 PetscCall(!cholmod_X_free_factor(&chol->factor,chol->common)); 273 } 274 if (chol->common->itype == CHOLMOD_INT) { 275 PetscCall(!cholmod_finish(chol->common)); 276 } else { 277 PetscCall(!cholmod_l_finish(chol->common)); 278 } 279 PetscCall(PetscFree(chol->common)); 280 PetscCall(PetscFree(chol->matrix)); 281 PetscCall(PetscObjectComposeFunction((PetscObject)F,"MatFactorGetSolverType_C",NULL)); 282 PetscCall(PetscObjectComposeFunction((PetscObject)F,"MatQRFactorSymbolic_C",NULL)); 283 PetscCall(PetscObjectComposeFunction((PetscObject)F,"MatQRFactorNumeric_C",NULL)); 284 PetscCall(PetscFree(F->data)); 285 PetscFunctionReturn(0); 286 } 287 288 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 289 static PetscErrorCode MatMatSolve_CHOLMOD(Mat,Mat,Mat); 290 291 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 292 293 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 294 { 295 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 296 const cholmod_common *c = chol->common; 297 PetscInt i; 298 299 PetscFunctionBegin; 300 if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(0); 301 PetscCall(PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n")); 302 PetscCall(PetscViewerASCIIPushTab(viewer)); 303 PetscCall(PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE")); 304 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound)); 305 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0)); 306 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1)); 307 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2)); 308 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank)); 309 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch)); 310 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal)); 311 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis)); 312 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super)); 313 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll)); 314 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack)); 315 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic)); 316 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol)); 317 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2])); 318 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2])); 319 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper)); 320 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print)); 321 for (i=0; i<c->nmethods; i++) { 322 PetscCall(PetscViewerASCIIPrintf(viewer,"Ordering method %" PetscInt_FMT "%s:\n",i,i==c->selected ? " [SELECTED]" : "")); 323 PetscCall(PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", 324 c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2)); 325 } 326 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder)); 327 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis)); 328 /* Statistics */ 329 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl)); 330 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz)); 331 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz)); 332 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl)); 333 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count)); 334 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage)); 335 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse)); 336 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col)); 337 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor)); 338 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit)); 339 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl)); 340 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl)); 341 #if defined(PETSC_USE_SUITESPARSE_GPU) 342 PetscCall(PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU)); 343 #endif 344 PetscCall(PetscViewerASCIIPopTab(viewer)); 345 PetscFunctionReturn(0); 346 } 347 348 PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 349 { 350 PetscBool iascii; 351 PetscViewerFormat format; 352 353 PetscFunctionBegin; 354 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 355 if (iascii) { 356 PetscCall(PetscViewerGetFormat(viewer,&format)); 357 if (format == PETSC_VIEWER_ASCII_INFO) { 358 PetscCall(MatView_Info_CHOLMOD(F,viewer)); 359 } 360 } 361 PetscFunctionReturn(0); 362 } 363 364 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 365 { 366 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 367 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 368 369 PetscFunctionBegin; 370 static_F = F; 371 PetscCall(VecWrapCholmod(B,GET_ARRAY_READ,&cholB)); 372 PetscCall(VecWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 373 X_handle = &cholX; 374 PetscCall(!cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common)); 375 PetscCall(!cholmod_X_free_dense(&Y_handle,chol->common)); 376 PetscCall(!cholmod_X_free_dense(&E_handle,chol->common)); 377 PetscCall(VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB)); 378 PetscCall(VecUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 379 PetscCall(PetscLogFlops(4.0*chol->common->lnz)); 380 PetscFunctionReturn(0); 381 } 382 383 static PetscErrorCode MatMatSolve_CHOLMOD(Mat F,Mat B,Mat X) 384 { 385 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 386 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 387 388 PetscFunctionBegin; 389 static_F = F; 390 PetscCall(MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB)); 391 PetscCall(MatDenseWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 392 X_handle = &cholX; 393 PetscCall(!cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common)); 394 PetscCall(!cholmod_X_free_dense(&Y_handle,chol->common)); 395 PetscCall(!cholmod_X_free_dense(&E_handle,chol->common)); 396 PetscCall(MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB)); 397 PetscCall(MatDenseUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX)); 398 PetscCall(PetscLogFlops(4.0*B->cmap->n*chol->common->lnz)); 399 PetscFunctionReturn(0); 400 } 401 402 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 403 { 404 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 405 cholmod_sparse cholA; 406 PetscBool aijalloc,valloc; 407 int err; 408 409 PetscFunctionBegin; 410 PetscCall((*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc)); 411 static_F = F; 412 err = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 413 PetscCheck(!err,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 414 PetscCheck(chol->common->status != CHOLMOD_NOT_POSDEF,PetscObjectComm((PetscObject)F),PETSC_ERR_MAT_CH_ZRPVT,"CHOLMOD detected that the matrix is not positive definite, failure at column %u",(unsigned)chol->factor->minor); 415 416 PetscCall(PetscLogFlops(chol->common->fl)); 417 if (aijalloc) PetscCall(PetscFree2(cholA.p,cholA.i)); 418 if (valloc) PetscCall(PetscFree(cholA.x)); 419 #if defined(PETSC_USE_SUITESPARSE_GPU) 420 PetscCall(PetscLogGpuTimeAdd(chol->common->CHOLMOD_GPU_GEMM_TIME + chol->common->CHOLMOD_GPU_SYRK_TIME + chol->common->CHOLMOD_GPU_TRSM_TIME + chol->common->CHOLMOD_GPU_POTRF_TIME)); 421 #endif 422 423 F->ops->solve = MatSolve_CHOLMOD; 424 F->ops->solvetranspose = MatSolve_CHOLMOD; 425 F->ops->matsolve = MatMatSolve_CHOLMOD; 426 F->ops->matsolvetranspose = MatMatSolve_CHOLMOD; 427 PetscFunctionReturn(0); 428 } 429 430 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 431 { 432 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 433 int err; 434 cholmod_sparse cholA; 435 PetscBool aijalloc,valloc; 436 PetscInt *fset = 0; 437 size_t fsize = 0; 438 439 PetscFunctionBegin; 440 PetscCall((*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc,&valloc)); 441 static_F = F; 442 if (chol->factor) { 443 err = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 444 PetscCheck(!err,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 445 } else if (perm) { 446 const PetscInt *ip; 447 PetscCall(ISGetIndices(perm,&ip)); 448 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 449 PetscCheck(chol->factor,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed using PETSc ordering with status %d",chol->common->status); 450 PetscCall(ISRestoreIndices(perm,&ip)); 451 } else { 452 chol->factor = cholmod_X_analyze(&cholA,chol->common); 453 PetscCheck(chol->factor,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed using internal ordering with status %d",chol->common->status); 454 } 455 456 if (aijalloc) PetscCall(PetscFree2(cholA.p,cholA.i)); 457 if (valloc) PetscCall(PetscFree(cholA.x)); 458 459 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 460 PetscFunctionReturn(0); 461 } 462 463 static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type) 464 { 465 PetscFunctionBegin; 466 *type = MATSOLVERCHOLMOD; 467 PetscFunctionReturn(0); 468 } 469 470 PETSC_INTERN PetscErrorCode MatGetInfo_CHOLMOD(Mat F,MatInfoType flag,MatInfo *info) 471 { 472 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 473 474 PetscFunctionBegin; 475 info->block_size = 1.0; 476 info->nz_allocated = chol->common->lnz; 477 info->nz_used = chol->common->lnz; 478 info->nz_unneeded = 0.0; 479 info->assemblies = 0.0; 480 info->mallocs = 0.0; 481 info->memory = chol->common->memory_inuse; 482 info->fill_ratio_given = 0; 483 info->fill_ratio_needed = 0; 484 info->factor_mallocs = chol->common->malloc_count; 485 PetscFunctionReturn(0); 486 } 487 488 /*MC 489 MATSOLVERCHOLMOD 490 491 A matrix type providing direct solvers (Cholesky) for sequential matrices 492 via the external package CHOLMOD. 493 494 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 495 496 Use -pc_type cholesky -pc_factor_mat_solver_type cholmod to use this direct solver 497 498 Consult CHOLMOD documentation for more information about the Common parameters 499 which correspond to the options database keys below. 500 501 Options Database Keys: 502 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 503 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 504 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 505 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 506 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 507 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 508 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 509 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 510 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 511 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 512 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 513 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 514 . -mat_cholmod_print <3> - Verbosity level (None) 515 - -mat_ordering_type internal - Use the ordering provided by Cholmod 516 517 Level: beginner 518 519 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 520 521 .seealso: `PCCHOLESKY`, `PCFactorSetMatSolverType()`, `MatSolverType` 522 M*/ 523 524 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 525 { 526 Mat B; 527 Mat_CHOLMOD *chol; 528 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 529 const char *prefix; 530 531 PetscFunctionBegin; 532 PetscCall(MatGetBlockSize(A,&bs)); 533 PetscCheck(bs == 1,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %" PetscInt_FMT,bs); 534 #if defined(PETSC_USE_COMPLEX) 535 PetscCheck(A->hermitian,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for Hermitian matrices"); 536 #endif 537 /* Create the factorization matrix F */ 538 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 539 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n)); 540 PetscCall(PetscStrallocpy("cholmod",&((PetscObject)B)->type_name)); 541 PetscCall(MatGetOptionsPrefix(A,&prefix)); 542 PetscCall(MatSetOptionsPrefix(B,prefix)); 543 PetscCall(MatSetUp(B)); 544 PetscCall(PetscNewLog(B,&chol)); 545 546 chol->Wrap = MatWrapCholmod_seqsbaij; 547 B->data = chol; 548 549 B->ops->getinfo = MatGetInfo_CHOLMOD; 550 B->ops->view = MatView_CHOLMOD; 551 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 552 B->ops->destroy = MatDestroy_CHOLMOD; 553 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod)); 554 B->factortype = MAT_FACTOR_CHOLESKY; 555 B->assembled = PETSC_TRUE; 556 B->preallocated = PETSC_TRUE; 557 558 PetscCall(CholmodStart(B)); 559 560 PetscCall(PetscFree(B->solvertype)); 561 PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype)); 562 B->canuseordering = PETSC_TRUE; 563 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_CHOLESKY])); 564 *F = B; 565 PetscFunctionReturn(0); 566 } 567