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 PetscErrorCode ierr; 37 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 38 cholmod_common *c; 39 PetscBool flg; 40 41 PetscFunctionBegin; 42 if (chol->common) PetscFunctionReturn(0); 43 PetscCall(PetscMalloc1(1,&chol->common)); 44 PetscCall(!cholmod_X_start(chol->common)); 45 46 c = chol->common; 47 c->error_handler = CholmodErrorHandler; 48 49 #define CHOLMOD_OPTION_DOUBLE(name,help) do { \ 50 PetscReal tmp = (PetscReal)c->name; \ 51 PetscCall(PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 52 c->name = (double)tmp; \ 53 } while (0) 54 55 #define CHOLMOD_OPTION_INT(name,help) do { \ 56 PetscInt tmp = (PetscInt)c->name; \ 57 PetscCall(PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 58 c->name = (int)tmp; \ 59 } while (0) 60 61 #define CHOLMOD_OPTION_SIZE_T(name,help) do { \ 62 PetscReal tmp = (PetscInt)c->name; \ 63 PetscCall(PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 64 PetscCheckFalse(tmp < 0,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"value must be positive"); \ 65 c->name = (size_t)tmp; \ 66 } while (0) 67 68 #define CHOLMOD_OPTION_BOOL(name,help) do { \ 69 PetscBool tmp = (PetscBool) !!c->name; \ 70 PetscCall(PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL)); \ 71 c->name = (int)tmp; \ 72 } while (0) 73 74 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat");PetscCall(ierr); 75 CHOLMOD_OPTION_INT(nmethods,"Number of different ordering methods to try"); 76 77 #if defined(PETSC_USE_SUITESPARSE_GPU) 78 c->useGPU = 1; 79 CHOLMOD_OPTION_INT(useGPU,"Use GPU for BLAS 1, otherwise 0"); 80 CHOLMOD_OPTION_SIZE_T(maxGpuMemBytes,"Maximum memory to allocate on the GPU"); 81 CHOLMOD_OPTION_DOUBLE(maxGpuMemFraction,"Fraction of available GPU memory to allocate"); 82 #endif 83 84 /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */ 85 chol->pack = (PetscBool)c->final_pack; 86 PetscCall(PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL)); 87 c->final_pack = (int)chol->pack; 88 89 CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D"); 90 CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified"); 91 CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified"); 92 CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified"); 93 CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]"); 94 { 95 static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0}; 96 PetscCall(PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL)); 97 } 98 if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization"); 99 CHOLMOD_OPTION_BOOL(final_asis,"Leave factors \"as is\""); 100 CHOLMOD_OPTION_BOOL(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)"); 101 if (!c->final_asis) { 102 CHOLMOD_OPTION_BOOL(final_super,"Leave supernodal factors instead of converting to simplicial"); 103 CHOLMOD_OPTION_BOOL(final_ll,"Turn LDL' factorization into LL'"); 104 CHOLMOD_OPTION_BOOL(final_monotonic,"Ensure columns are monotonic when done"); 105 CHOLMOD_OPTION_BOOL(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation"); 106 } 107 { 108 PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]}; 109 PetscInt n = 3; 110 PetscCall(PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg)); 111 PetscCheckFalse(flg && n != 3,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_zrelax"); 112 if (flg) while (n--) c->zrelax[n] = (double)tmp[n]; 113 } 114 { 115 PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]}; 116 PetscCall(PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg)); 117 PetscCheckFalse(flg && n != 3,PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_nrelax"); 118 if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n]; 119 } 120 CHOLMOD_OPTION_BOOL(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]"); 121 CHOLMOD_OPTION_BOOL(default_nesdis,"Use NESDIS instead of METIS for nested dissection"); 122 CHOLMOD_OPTION_INT(print,"Verbosity level"); 123 ierr = PetscOptionsEnd();PetscCall(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc,PetscBool *valloc) 128 { 129 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)A->data; 130 PetscBool vallocin = PETSC_FALSE; 131 132 PetscFunctionBegin; 133 PetscCall(PetscMemzero(C,sizeof(*C))); 134 /* 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 */ 135 C->nrow = (size_t)A->cmap->n; 136 C->ncol = (size_t)A->rmap->n; 137 C->nzmax = (size_t)sbaij->maxnz; 138 C->p = sbaij->i; 139 C->i = sbaij->j; 140 if (values) { 141 #if defined(PETSC_USE_COMPLEX) 142 /* we need to pass CHOLMOD the conjugate matrix */ 143 PetscScalar *v; 144 PetscInt i; 145 146 PetscCall(PetscMalloc1(sbaij->maxnz,&v)); 147 for (i = 0; i < sbaij->maxnz; i++) v[i] = PetscConj(sbaij->a[i]); 148 C->x = v; 149 vallocin = PETSC_TRUE; 150 #else 151 C->x = sbaij->a; 152 #endif 153 } 154 C->stype = -1; 155 C->itype = CHOLMOD_INT_TYPE; 156 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 157 C->dtype = CHOLMOD_DOUBLE; 158 C->sorted = 1; 159 C->packed = 1; 160 *aijalloc = PETSC_FALSE; 161 *valloc = vallocin; 162 PetscFunctionReturn(0); 163 } 164 165 #define GET_ARRAY_READ 0 166 #define GET_ARRAY_WRITE 1 167 168 PetscErrorCode VecWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y) 169 { 170 PetscScalar *x; 171 PetscInt n; 172 173 PetscFunctionBegin; 174 PetscCall(PetscMemzero(Y,sizeof(*Y))); 175 switch (rw) { 176 case GET_ARRAY_READ: 177 PetscCall(VecGetArrayRead(X,(const PetscScalar**)&x)); 178 break; 179 case GET_ARRAY_WRITE: 180 PetscCall(VecGetArrayWrite(X,&x)); 181 break; 182 default: 183 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 184 break; 185 } 186 PetscCall(VecGetSize(X,&n)); 187 188 Y->x = x; 189 Y->nrow = n; 190 Y->ncol = 1; 191 Y->nzmax = n; 192 Y->d = n; 193 Y->xtype = CHOLMOD_SCALAR_TYPE; 194 Y->dtype = CHOLMOD_DOUBLE; 195 PetscFunctionReturn(0); 196 } 197 198 PetscErrorCode VecUnWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y) 199 { 200 PetscFunctionBegin; 201 switch (rw) { 202 case GET_ARRAY_READ: 203 PetscCall(VecRestoreArrayRead(X,(const PetscScalar**)&Y->x)); 204 break; 205 case GET_ARRAY_WRITE: 206 PetscCall(VecRestoreArrayWrite(X,(PetscScalar**)&Y->x)); 207 break; 208 default: 209 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 210 break; 211 } 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode MatDenseWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 216 { 217 PetscScalar *x; 218 PetscInt m,n,lda; 219 220 PetscFunctionBegin; 221 PetscCall(PetscMemzero(Y,sizeof(*Y))); 222 switch (rw) { 223 case GET_ARRAY_READ: 224 PetscCall(MatDenseGetArrayRead(X,(const PetscScalar**)&x)); 225 break; 226 case GET_ARRAY_WRITE: 227 PetscCall(MatDenseGetArrayWrite(X,&x)); 228 break; 229 default: 230 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 231 break; 232 } 233 PetscCall(MatDenseGetLDA(X,&lda)); 234 PetscCall(MatGetLocalSize(X,&m,&n)); 235 236 Y->x = x; 237 Y->nrow = m; 238 Y->ncol = n; 239 Y->nzmax = lda*n; 240 Y->d = lda; 241 Y->xtype = CHOLMOD_SCALAR_TYPE; 242 Y->dtype = CHOLMOD_DOUBLE; 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode MatDenseUnWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 247 { 248 PetscFunctionBegin; 249 switch (rw) { 250 case GET_ARRAY_READ: 251 PetscCall(MatDenseRestoreArrayRead(X,(const PetscScalar**)&Y->x)); 252 break; 253 case GET_ARRAY_WRITE: 254 /* we don't have MatDenseRestoreArrayWrite */ 255 PetscCall(MatDenseRestoreArray(X,(PetscScalar**)&Y->x)); 256 break; 257 default: 258 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw); 259 break; 260 } 261 PetscFunctionReturn(0); 262 } 263 264 PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F) 265 { 266 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 267 268 PetscFunctionBegin; 269 if (chol->spqrfact) { 270 PetscCall(!SuiteSparseQR_C_free(&chol->spqrfact, chol->common)); 271 } 272 if (chol->factor) { 273 PetscCall(!cholmod_X_free_factor(&chol->factor,chol->common)); 274 } 275 if (chol->common->itype == CHOLMOD_INT) { 276 PetscCall(!cholmod_finish(chol->common)); 277 } else { 278 PetscCall(!cholmod_l_finish(chol->common)); 279 } 280 PetscCall(PetscFree(chol->common)); 281 PetscCall(PetscFree(chol->matrix)); 282 PetscCall(PetscObjectComposeFunction((PetscObject)F,"MatFactorGetSolverType_C",NULL)); 283 PetscCall(PetscFree(F->data)); 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 288 static PetscErrorCode MatMatSolve_CHOLMOD(Mat,Mat,Mat); 289 290 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 291 292 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 293 { 294 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 295 const cholmod_common *c = chol->common; 296 PetscErrorCode ierr; 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 ierr = 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);PetscCall(ierr); 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 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 PetscCall((*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc)); 411 static_F = F; 412 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 413 PetscCheck(!ierr,PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 414 PetscCheckFalse(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 PetscErrorCode ierr; 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 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 444 PetscCheck(!ierr,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 PetscCheckFalse(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