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