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 static 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 static 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 static 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 /* we don't have MatDenseGetArrayWrite */ 235 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 236 break; 237 default: 238 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D not handled",rw); 239 break; 240 } 241 ierr = MatDenseGetLDA(X,&lda);CHKERRQ(ierr); 242 ierr = MatGetLocalSize(X,&m,&n);CHKERRQ(ierr); 243 244 Y->x = x; 245 Y->nrow = m; 246 Y->ncol = n; 247 Y->nzmax = lda*n; 248 Y->d = lda; 249 Y->xtype = CHOLMOD_SCALAR_TYPE; 250 Y->dtype = CHOLMOD_DOUBLE; 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode MatDenseUnWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y) 255 { 256 PetscErrorCode ierr; 257 258 PetscFunctionBegin; 259 switch (rw) { 260 case GET_ARRAY_READ: 261 ierr = MatDenseRestoreArrayRead(X,(const PetscScalar**)&Y->x);CHKERRQ(ierr); 262 break; 263 case GET_ARRAY_WRITE: 264 /* we don't have MatDenseRestoreArrayWrite */ 265 ierr = MatDenseRestoreArray(X,(PetscScalar**)&Y->x);CHKERRQ(ierr); 266 break; 267 default: 268 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D not handled",rw); 269 break; 270 } 271 PetscFunctionReturn(0); 272 } 273 274 PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F) 275 { 276 PetscErrorCode ierr; 277 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 278 279 PetscFunctionBegin; 280 ierr = !cholmod_X_free_factor(&chol->factor,chol->common);CHKERRQ(ierr); 281 ierr = !cholmod_X_finish(chol->common);CHKERRQ(ierr); 282 ierr = PetscFree(chol->common);CHKERRQ(ierr); 283 ierr = PetscFree(chol->matrix);CHKERRQ(ierr); 284 ierr = PetscObjectComposeFunction((PetscObject)F,"MatFactorGetSolverType_C",NULL);CHKERRQ(ierr); 285 ierr = PetscFree(F->data);CHKERRQ(ierr); 286 PetscFunctionReturn(0); 287 } 288 289 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 290 static PetscErrorCode MatMatSolve_CHOLMOD(Mat,Mat,Mat); 291 292 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 293 294 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 295 { 296 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 297 const cholmod_common *c = chol->common; 298 PetscErrorCode ierr; 299 PetscInt i; 300 301 PetscFunctionBegin; 302 if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(0); 303 ierr = PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n");CHKERRQ(ierr); 304 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 305 ierr = PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE");CHKERRQ(ierr); 306 ierr = PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound);CHKERRQ(ierr); 307 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0);CHKERRQ(ierr); 308 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1);CHKERRQ(ierr); 309 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2);CHKERRQ(ierr); 310 ierr = PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank);CHKERRQ(ierr); 311 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch);CHKERRQ(ierr); 312 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal);CHKERRQ(ierr); 313 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis);CHKERRQ(ierr); 314 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super);CHKERRQ(ierr); 315 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll);CHKERRQ(ierr); 316 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack);CHKERRQ(ierr); 317 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic);CHKERRQ(ierr); 318 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2]);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2]);CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print);CHKERRQ(ierr); 323 for (i=0; i<c->nmethods; i++) { 324 ierr = PetscViewerASCIIPrintf(viewer,"Ordering method %D%s:\n",i,i==c->selected ? " [SELECTED]" : "");CHKERRQ(ierr); 325 ierr = PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", 326 c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2);CHKERRQ(ierr); 327 } 328 ierr = PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder);CHKERRQ(ierr); 329 ierr = PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis);CHKERRQ(ierr); 330 /* Statistics */ 331 ierr = PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl);CHKERRQ(ierr); 332 ierr = PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz);CHKERRQ(ierr); 333 ierr = PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz);CHKERRQ(ierr); 334 ierr = PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl);CHKERRQ(ierr); 335 ierr = PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count);CHKERRQ(ierr);CHKERRQ(ierr); 336 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage);CHKERRQ(ierr);CHKERRQ(ierr); 337 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse);CHKERRQ(ierr);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col);CHKERRQ(ierr);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor);CHKERRQ(ierr);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit);CHKERRQ(ierr);CHKERRQ(ierr); 341 ierr = PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl);CHKERRQ(ierr);CHKERRQ(ierr); 342 ierr = PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl);CHKERRQ(ierr);CHKERRQ(ierr); 343 #if defined(PETSC_USE_SUITESPARSE_GPU) 344 ierr = PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU);CHKERRQ(ierr);CHKERRQ(ierr); 345 #endif 346 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 350 PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 351 { 352 PetscErrorCode ierr; 353 PetscBool iascii; 354 PetscViewerFormat format; 355 356 PetscFunctionBegin; 357 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 358 if (iascii) { 359 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 360 if (format == PETSC_VIEWER_ASCII_INFO) { 361 ierr = MatView_Info_CHOLMOD(F,viewer);CHKERRQ(ierr); 362 } 363 } 364 PetscFunctionReturn(0); 365 } 366 367 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 368 { 369 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 370 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 371 PetscErrorCode ierr; 372 373 PetscFunctionBegin; 374 static_F = F; 375 ierr = VecWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 376 ierr = VecWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 377 X_handle = &cholX; 378 ierr = !cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common);CHKERRQ(ierr); 379 ierr = !cholmod_X_free_dense(&Y_handle,chol->common);CHKERRQ(ierr); 380 ierr = !cholmod_X_free_dense(&E_handle,chol->common);CHKERRQ(ierr); 381 ierr = VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 382 ierr = VecUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 383 PetscFunctionReturn(0); 384 } 385 386 static PetscErrorCode MatMatSolve_CHOLMOD(Mat F,Mat B,Mat X) 387 { 388 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 389 cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL; 390 PetscErrorCode ierr; 391 392 PetscFunctionBegin; 393 static_F = F; 394 ierr = MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 395 ierr = MatDenseWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 396 X_handle = &cholX; 397 ierr = !cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common);CHKERRQ(ierr); 398 ierr = !cholmod_X_free_dense(&Y_handle,chol->common);CHKERRQ(ierr); 399 ierr = !cholmod_X_free_dense(&E_handle,chol->common);CHKERRQ(ierr); 400 ierr = MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);CHKERRQ(ierr); 401 ierr = MatDenseUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 406 { 407 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 408 cholmod_sparse cholA; 409 PetscBool aijalloc,valloc; 410 PetscErrorCode ierr; 411 412 PetscFunctionBegin; 413 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);CHKERRQ(ierr); 414 static_F = F; 415 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 416 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 417 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); 418 419 if (aijalloc) {ierr = PetscFree2(cholA.p,cholA.i);CHKERRQ(ierr);} 420 if (valloc) {ierr = PetscFree(cholA.x);CHKERRQ(ierr);} 421 #if defined(PETSC_USE_SUITESPARSE_GPU) 422 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); 423 #endif 424 425 F->ops->solve = MatSolve_CHOLMOD; 426 F->ops->solvetranspose = MatSolve_CHOLMOD; 427 F->ops->matsolve = MatMatSolve_CHOLMOD; 428 F->ops->matsolvetranspose = MatMatSolve_CHOLMOD; 429 PetscFunctionReturn(0); 430 } 431 432 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 433 { 434 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 435 PetscErrorCode ierr; 436 cholmod_sparse cholA; 437 PetscBool aijalloc,valloc; 438 PetscInt *fset = 0; 439 size_t fsize = 0; 440 441 PetscFunctionBegin; 442 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc,&valloc);CHKERRQ(ierr); 443 static_F = F; 444 if (chol->factor) { 445 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 446 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 447 } else if (perm) { 448 const PetscInt *ip; 449 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 450 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 451 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 452 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 453 } else { 454 chol->factor = cholmod_X_analyze(&cholA,chol->common); 455 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 456 } 457 458 if (aijalloc) {ierr = PetscFree2(cholA.p,cholA.i);CHKERRQ(ierr);} 459 if (valloc) {ierr = PetscFree(cholA.x);CHKERRQ(ierr);} 460 461 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 462 PetscFunctionReturn(0); 463 } 464 465 static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type) 466 { 467 PetscFunctionBegin; 468 *type = MATSOLVERCHOLMOD; 469 PetscFunctionReturn(0); 470 } 471 472 PETSC_INTERN PetscErrorCode MatGetInfo_CHOLMOD(Mat F,MatInfoType flag,MatInfo *info) 473 { 474 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 475 476 PetscFunctionBegin; 477 info->block_size = 1.0; 478 info->nz_allocated = chol->common->lnz; 479 info->nz_used = chol->common->lnz; 480 info->nz_unneeded = 0.0; 481 info->assemblies = 0.0; 482 info->mallocs = 0.0; 483 info->memory = chol->common->memory_inuse; 484 info->fill_ratio_given = 0; 485 info->fill_ratio_needed = 0; 486 info->factor_mallocs = chol->common->malloc_count; 487 PetscFunctionReturn(0); 488 } 489 490 /*MC 491 MATSOLVERCHOLMOD = "cholmod" - 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 516 Level: beginner 517 518 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 519 520 .seealso: PCCHOLESKY, PCFactorSetMatSolverType(), MatSolverType 521 M*/ 522 523 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 524 { 525 Mat B; 526 Mat_CHOLMOD *chol; 527 PetscErrorCode ierr; 528 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 529 const char *prefix; 530 531 PetscFunctionBegin; 532 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 533 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 534 #if defined(PETSC_USE_COMPLEX) 535 if (!A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for hermitian matrices"); 536 #endif 537 /* Create the factorization matrix F */ 538 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 539 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 540 ierr = PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);CHKERRQ(ierr); 541 ierr = MatGetOptionsPrefix(A,&prefix);CHKERRQ(ierr); 542 ierr = MatSetOptionsPrefix(B,prefix);CHKERRQ(ierr); 543 ierr = MatSetUp(B);CHKERRQ(ierr); 544 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 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 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod);CHKERRQ(ierr); 554 B->factortype = MAT_FACTOR_CHOLESKY; 555 B->assembled = PETSC_TRUE; 556 B->preallocated = PETSC_TRUE; 557 558 ierr = CholmodStart(B);CHKERRQ(ierr); 559 560 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 561 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 562 563 *F = B; 564 PetscFunctionReturn(0); 565 } 566