1 2 #include <petscsys.h> 3 #include <../src/mat/impls/aij/seq/aij.h> 4 #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h> 5 6 EXTERN_C_BEGIN 7 #include <SuiteSparseQR_C.h> 8 EXTERN_C_END 9 10 static PetscErrorCode MatWrapCholmod_SPQR_seqaij(Mat A, PetscBool values, cholmod_sparse *C, PetscBool *aijalloc, PetscBool *valloc) { 11 Mat_SeqAIJ *aij; 12 Mat AT; 13 const PetscScalar *aa; 14 PetscScalar *ca; 15 const PetscInt *ai, *aj; 16 PetscInt n = A->cmap->n, i, j, k, nz; 17 SuiteSparse_long *ci, *cj; /* SuiteSparse_long is the only choice for SPQR */ 18 PetscBool vain = PETSC_FALSE, flg; 19 20 PetscFunctionBegin; 21 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg)); 22 if (flg) { 23 PetscCall(MatNormalHermitianGetMat(A, &A)); 24 } else if (!PetscDefined(USE_COMPLEX)) { 25 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg)); 26 if (flg) { PetscCall(MatNormalGetMat(A, &A)); } 27 } 28 /* cholmod_sparse is compressed sparse column */ 29 PetscCall(MatIsSymmetric(A, 0.0, &flg)); 30 if (flg) { 31 PetscCall(PetscObjectReference((PetscObject)A)); 32 AT = A; 33 } else { 34 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT)); 35 } 36 aij = (Mat_SeqAIJ *)AT->data; 37 ai = aij->j; 38 aj = aij->i; 39 for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j]; 40 PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci)); 41 if (values) { 42 vain = PETSC_TRUE; 43 PetscCall(PetscMalloc1(nz, &ca)); 44 PetscCall(MatSeqAIJGetArrayRead(AT, &aa)); 45 } 46 for (j = 0, k = 0; j < n; j++) { 47 cj[j] = k; 48 for (i = aj[j]; i < aj[j + 1]; i++, k++) { 49 ci[k] = ai[i]; 50 if (values) ca[k] = aa[i]; 51 } 52 } 53 cj[j] = k; 54 *aijalloc = PETSC_TRUE; 55 *valloc = vain; 56 if (values) { PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa)); } 57 58 PetscCall(PetscMemzero(C, sizeof(*C))); 59 60 C->nrow = (size_t)AT->cmap->n; 61 C->ncol = (size_t)AT->rmap->n; 62 C->nzmax = (size_t)nz; 63 C->p = cj; 64 C->i = ci; 65 C->x = values ? ca : 0; 66 C->stype = 0; 67 C->itype = CHOLMOD_LONG; 68 C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN; 69 C->dtype = CHOLMOD_DOUBLE; 70 C->sorted = 1; 71 C->packed = 1; 72 73 PetscCall(MatDestroy(&AT)); 74 PetscFunctionReturn(0); 75 } 76 77 static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type) { 78 PetscFunctionBegin; 79 *type = MATSOLVERSPQR; 80 PetscFunctionReturn(0); 81 } 82 83 #define GET_ARRAY_READ 0 84 #define GET_ARRAY_WRITE 1 85 86 static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle) { 87 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 88 cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL; 89 90 PetscFunctionBegin; 91 if (!chol->normal) { 92 QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common); 93 PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed"); 94 Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common); 95 PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 96 } else { 97 Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common); 98 PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 99 Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common); 100 PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 101 PetscCall(!cholmod_l_free_dense(&Z_handle, chol->common)); 102 } 103 *_Y_handle = Y_handle; 104 PetscCall(!cholmod_l_free_dense(&QTB_handle, chol->common)); 105 PetscFunctionReturn(0); 106 } 107 108 static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X) { 109 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 110 cholmod_dense cholB, *Y_handle = NULL; 111 PetscInt n; 112 PetscScalar *v; 113 114 PetscFunctionBegin; 115 PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB)); 116 PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle)); 117 PetscCall(VecGetLocalSize(X, &n)); 118 PetscCall(VecGetArrayWrite(X, &v)); 119 PetscCall(PetscArraycpy(v, (PetscScalar *)(Y_handle->x), n)); 120 PetscCall(VecRestoreArrayWrite(X, &v)); 121 PetscCall(!cholmod_l_free_dense(&Y_handle, chol->common)); 122 PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB)); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X) { 127 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 128 cholmod_dense cholB, *Y_handle = NULL; 129 PetscScalar *v; 130 PetscInt lda; 131 132 PetscFunctionBegin; 133 PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB)); 134 PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle)); 135 PetscCall(MatDenseGetArrayWrite(X, &v)); 136 PetscCall(MatDenseGetLDA(X, &lda)); 137 if ((size_t)lda == Y_handle->d) { 138 PetscCall(PetscArraycpy(v, (PetscScalar *)(Y_handle->x), lda * Y_handle->ncol)); 139 } else { 140 for (size_t j = 0; j < Y_handle->ncol; j++) { PetscCall(PetscArraycpy(&v[j * lda], &(((PetscScalar *)Y_handle->x)[j * Y_handle->d]), Y_handle->nrow)); } 141 } 142 PetscCall(MatDenseRestoreArrayWrite(X, &v)); 143 PetscCall(!cholmod_l_free_dense(&Y_handle, chol->common)); 144 PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB)); 145 PetscFunctionReturn(0); 146 } 147 148 static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle) { 149 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 150 cholmod_dense *Y_handle = NULL, *RTB_handle = NULL; 151 152 PetscFunctionBegin; 153 RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common); 154 PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed"); 155 Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common); 156 PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed"); 157 *_Y_handle = Y_handle; 158 PetscCall(!cholmod_l_free_dense(&RTB_handle, chol->common)); 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X) { 163 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 164 cholmod_dense cholB, *Y_handle = NULL; 165 PetscInt n; 166 PetscScalar *v; 167 168 PetscFunctionBegin; 169 PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB)); 170 PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle)); 171 PetscCall(VecGetLocalSize(X, &n)); 172 PetscCall(VecGetArrayWrite(X, &v)); 173 PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n)); 174 PetscCall(VecRestoreArrayWrite(X, &v)); 175 PetscCall(!cholmod_l_free_dense(&Y_handle, chol->common)); 176 PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB)); 177 PetscFunctionReturn(0); 178 } 179 180 static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X) { 181 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 182 cholmod_dense cholB, *Y_handle = NULL; 183 PetscScalar *v; 184 PetscInt lda; 185 186 PetscFunctionBegin; 187 PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB)); 188 PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle)); 189 PetscCall(MatDenseGetArrayWrite(X, &v)); 190 PetscCall(MatDenseGetLDA(X, &lda)); 191 if ((size_t)lda == Y_handle->d) { 192 PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol)); 193 } else { 194 for (size_t j = 0; j < Y_handle->ncol; j++) { PetscCall(PetscArraycpy(&v[j * lda], &(((PetscScalar *)Y_handle->x)[j * Y_handle->d]), Y_handle->nrow)); } 195 } 196 PetscCall(MatDenseRestoreArrayWrite(X, &v)); 197 PetscCall(!cholmod_l_free_dense(&Y_handle, chol->common)); 198 PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB)); 199 PetscFunctionReturn(0); 200 } 201 202 static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info) { 203 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 204 cholmod_sparse cholA; 205 PetscBool aijalloc, valloc; 206 int err; 207 208 PetscFunctionBegin; 209 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal)); 210 if (!chol->normal && !PetscDefined(USE_COMPLEX)) { PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal)); } 211 PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc)); 212 err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common); 213 PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status); 214 215 if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i)); 216 if (valloc) PetscCall(PetscFree(cholA.x)); 217 218 F->ops->solve = MatSolve_SPQR; 219 F->ops->matsolve = MatMatSolve_SPQR; 220 if (chol->normal) { 221 F->ops->solvetranspose = MatSolve_SPQR; 222 F->ops->matsolvetranspose = MatMatSolve_SPQR; 223 } else if (A->cmap->n == A->rmap->n) { 224 F->ops->solvetranspose = MatSolveTranspose_SPQR; 225 F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR; 226 } 227 PetscFunctionReturn(0); 228 } 229 230 PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info) { 231 Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data; 232 cholmod_sparse cholA; 233 PetscBool aijalloc, valloc; 234 235 PetscFunctionBegin; 236 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal)); 237 if (!chol->normal && !PetscDefined(USE_COMPLEX)) { PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal)); } 238 PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc)); 239 if (PetscDefined(USE_DEBUG)) { PetscCall(!cholmod_l_check_sparse(&cholA, chol->common)); } 240 if (chol->spqrfact) { PetscCall(!SuiteSparseQR_C_free(&chol->spqrfact, chol->common)); } 241 chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common); 242 PetscCheck(chol->spqrfact, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status); 243 244 if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i)); 245 if (valloc) PetscCall(PetscFree(cholA.x)); 246 247 PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR)); 248 PetscFunctionReturn(0); 249 } 250 251 /*MC 252 MATSOLVERSPQR 253 254 A matrix type providing direct solvers (QR factorizations) for sequential matrices 255 via the external package SPQR. 256 257 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 258 259 Consult SPQR documentation for more information about the Common parameters 260 which correspond to the options database keys below. 261 262 Level: beginner 263 264 Note: SPQR is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 265 266 .seealso: `PCQR`, `PCFactorSetMatSolverType()`, `MatSolverType` 267 M*/ 268 269 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A, MatFactorType ftype, Mat *F) { 270 Mat B; 271 Mat_CHOLMOD *chol; 272 PetscInt m = A->rmap->n, n = A->cmap->n; 273 const char *prefix; 274 275 PetscFunctionBegin; 276 /* Create the factorization matrix F */ 277 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 278 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 279 PetscCall(PetscStrallocpy("spqr", &((PetscObject)B)->type_name)); 280 PetscCall(MatGetOptionsPrefix(A, &prefix)); 281 PetscCall(MatSetOptionsPrefix(B, prefix)); 282 PetscCall(MatSetUp(B)); 283 PetscCall(PetscNewLog(B, &chol)); 284 285 chol->Wrap = MatWrapCholmod_SPQR_seqaij; 286 B->data = chol; 287 288 B->ops->getinfo = MatGetInfo_CHOLMOD; 289 B->ops->view = MatView_CHOLMOD; 290 B->ops->destroy = MatDestroy_CHOLMOD; 291 292 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_SPQR)); 293 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR)); 294 295 B->factortype = MAT_FACTOR_QR; 296 B->assembled = PETSC_TRUE; 297 B->preallocated = PETSC_TRUE; 298 299 PetscCall(PetscFree(B->solvertype)); 300 PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype)); 301 B->canuseordering = PETSC_FALSE; 302 PetscCall(CholmodStart(B)); 303 chol->common->itype = CHOLMOD_LONG; 304 *F = B; 305 PetscFunctionReturn(0); 306 } 307