Actual source code: aijspqr.c


  2: #include <petscsys.h>
  3: #include <../src/mat/impls/aij/seq/aij.h>
  4: #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h>

  6: EXTERN_C_BEGIN
  7: #include <SuiteSparseQR_C.h>
  8: EXTERN_C_END

 10: static PetscErrorCode MatWrapCholmod_SPQR_seqaij(Mat A, PetscBool values, cholmod_sparse *C, PetscBool *aijalloc, PetscBool *valloc)
 11: {
 12:   Mat_SeqAIJ        *aij;
 13:   Mat                AT;
 14:   const PetscScalar *aa;
 15:   PetscScalar       *ca;
 16:   const PetscInt    *ai, *aj;
 17:   PetscInt           n = A->cmap->n, i, j, k, nz;
 18:   SuiteSparse_long  *ci, *cj; /* SuiteSparse_long is the only choice for SPQR */
 19:   PetscBool          vain = PETSC_FALSE, flg;

 21:   PetscFunctionBegin;
 22:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg));
 23:   if (flg) {
 24:     PetscCall(MatNormalHermitianGetMat(A, &A));
 25:   } else if (!PetscDefined(USE_COMPLEX)) {
 26:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg));
 27:     if (flg) PetscCall(MatNormalGetMat(A, &A));
 28:   }
 29:   /* cholmod_sparse is compressed sparse column */
 30:   PetscCall(MatIsSymmetric(A, 0.0, &flg));
 31:   if (flg) {
 32:     PetscCall(PetscObjectReference((PetscObject)A));
 33:     AT = A;
 34:   } else {
 35:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
 36:   }
 37:   aij = (Mat_SeqAIJ *)AT->data;
 38:   ai  = aij->j;
 39:   aj  = aij->i;
 40:   for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
 41:   PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci));
 42:   if (values) {
 43:     vain = PETSC_TRUE;
 44:     PetscCall(PetscMalloc1(nz, &ca));
 45:     PetscCall(MatSeqAIJGetArrayRead(AT, &aa));
 46:   }
 47:   for (j = 0, k = 0; j < n; j++) {
 48:     cj[j] = k;
 49:     for (i = aj[j]; i < aj[j + 1]; i++, k++) {
 50:       ci[k] = ai[i];
 51:       if (values) ca[k] = aa[i];
 52:     }
 53:   }
 54:   cj[j]     = k;
 55:   *aijalloc = PETSC_TRUE;
 56:   *valloc   = vain;
 57:   if (values) PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa));

 59:   PetscCall(PetscMemzero(C, sizeof(*C)));

 61:   C->nrow   = (size_t)AT->cmap->n;
 62:   C->ncol   = (size_t)AT->rmap->n;
 63:   C->nzmax  = (size_t)nz;
 64:   C->p      = cj;
 65:   C->i      = ci;
 66:   C->x      = values ? ca : 0;
 67:   C->stype  = 0;
 68:   C->itype  = CHOLMOD_LONG;
 69:   C->xtype  = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
 70:   C->dtype  = CHOLMOD_DOUBLE;
 71:   C->sorted = 1;
 72:   C->packed = 1;

 74:   PetscCall(MatDestroy(&AT));
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
 79: {
 80:   PetscFunctionBegin;
 81:   *type = MATSOLVERSPQR;
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: #define GET_ARRAY_READ  0
 86: #define GET_ARRAY_WRITE 1

 88: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
 89: {
 90:   Mat_CHOLMOD   *chol     = (Mat_CHOLMOD *)F->data;
 91:   cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;

 93:   PetscFunctionBegin;
 94:   if (!chol->normal) {
 95:     QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
 96:     PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
 97:     Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
 98:     PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
 99:   } else {
100:     Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
101:     PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
102:     Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
103:     PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
104:     PetscCallExternal(!cholmod_l_free_dense, &Z_handle, chol->common);
105:   }
106:   *_Y_handle = Y_handle;
107:   PetscCallExternal(!cholmod_l_free_dense, &QTB_handle, chol->common);
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X)
112: {
113:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
114:   cholmod_dense cholB, *Y_handle = NULL;
115:   PetscInt      n;
116:   PetscScalar  *v;

118:   PetscFunctionBegin;
119:   PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
120:   PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
121:   PetscCall(VecGetLocalSize(X, &n));
122:   PetscCall(VecGetArrayWrite(X, &v));
123:   PetscCall(PetscArraycpy(v, (PetscScalar *)(Y_handle->x), n));
124:   PetscCall(VecRestoreArrayWrite(X, &v));
125:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
126:   PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
131: {
132:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
133:   cholmod_dense cholB, *Y_handle = NULL;
134:   PetscScalar  *v;
135:   PetscInt      lda;

137:   PetscFunctionBegin;
138:   PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
139:   PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
140:   PetscCall(MatDenseGetArrayWrite(X, &v));
141:   PetscCall(MatDenseGetLDA(X, &lda));
142:   if ((size_t)lda == Y_handle->d) {
143:     PetscCall(PetscArraycpy(v, (PetscScalar *)(Y_handle->x), lda * Y_handle->ncol));
144:   } else {
145:     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));
146:   }
147:   PetscCall(MatDenseRestoreArrayWrite(X, &v));
148:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
149:   PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
154: {
155:   Mat_CHOLMOD   *chol     = (Mat_CHOLMOD *)F->data;
156:   cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;

158:   PetscFunctionBegin;
159:   RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
160:   PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
161:   Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
162:   PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
163:   *_Y_handle = Y_handle;
164:   PetscCallExternal(!cholmod_l_free_dense, &RTB_handle, chol->common);
165:   PetscFunctionReturn(PETSC_SUCCESS);
166: }

168: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
169: {
170:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
171:   cholmod_dense cholB, *Y_handle = NULL;
172:   PetscInt      n;
173:   PetscScalar  *v;

175:   PetscFunctionBegin;
176:   PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
177:   PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
178:   PetscCall(VecGetLocalSize(X, &n));
179:   PetscCall(VecGetArrayWrite(X, &v));
180:   PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
181:   PetscCall(VecRestoreArrayWrite(X, &v));
182:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
183:   PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
188: {
189:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
190:   cholmod_dense cholB, *Y_handle = NULL;
191:   PetscScalar  *v;
192:   PetscInt      lda;

194:   PetscFunctionBegin;
195:   PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
196:   PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
197:   PetscCall(MatDenseGetArrayWrite(X, &v));
198:   PetscCall(MatDenseGetLDA(X, &lda));
199:   if ((size_t)lda == Y_handle->d) {
200:     PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
201:   } else {
202:     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));
203:   }
204:   PetscCall(MatDenseRestoreArrayWrite(X, &v));
205:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
206:   PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
211: {
212:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
213:   cholmod_sparse cholA;
214:   PetscBool      aijalloc, valloc;
215:   int            err;

217:   PetscFunctionBegin;
218:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
219:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
220:   PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
221:   err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);
222:   PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status);

224:   if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
225:   if (valloc) PetscCall(PetscFree(cholA.x));

227:   F->ops->solve    = MatSolve_SPQR;
228:   F->ops->matsolve = MatMatSolve_SPQR;
229:   if (chol->normal) {
230:     F->ops->solvetranspose    = MatSolve_SPQR;
231:     F->ops->matsolvetranspose = MatMatSolve_SPQR;
232:   } else if (A->cmap->n == A->rmap->n) {
233:     F->ops->solvetranspose    = MatSolveTranspose_SPQR;
234:     F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR;
235:   }
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
240: {
241:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
242:   cholmod_sparse cholA;
243:   PetscBool      aijalloc, valloc;

245:   PetscFunctionBegin;
246:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
247:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
248:   PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
249:   if (PetscDefined(USE_DEBUG)) PetscCallExternal(!cholmod_l_check_sparse, &cholA, chol->common);
250:   if (chol->spqrfact) PetscCallExternal(!SuiteSparseQR_C_free, &chol->spqrfact, chol->common);
251:   chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common);
252:   PetscCheck(chol->spqrfact, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status);

254:   if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
255:   if (valloc) PetscCall(PetscFree(cholA.x));

257:   PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR));
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: /*MC
262:   MATSOLVERSPQR

264:   A matrix type providing direct solvers, QR factorizations, for sequential matrices
265:   via the external package SPQR.

267:   Use `./configure --download-suitesparse` to install PETSc to use SPQR

269:   Consult SPQR documentation for more information about the common parameters
270:   which correspond to the options database keys below.

272:    Level: beginner

274:    Note:
275:    SPQR is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html

277: .seealso: [](ch_matrices), `Mat`, `PCQR`, `PCFactorSetMatSolverType()`, `MatSolverType`
278: M*/

280: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A, MatFactorType ftype, Mat *F)
281: {
282:   Mat          B;
283:   Mat_CHOLMOD *chol;
284:   PetscInt     m = A->rmap->n, n = A->cmap->n;
285:   const char  *prefix;

287:   PetscFunctionBegin;
288:   /* Create the factorization matrix F */
289:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
290:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
291:   PetscCall(PetscStrallocpy("spqr", &((PetscObject)B)->type_name));
292:   PetscCall(MatGetOptionsPrefix(A, &prefix));
293:   PetscCall(MatSetOptionsPrefix(B, prefix));
294:   PetscCall(MatSetUp(B));
295:   PetscCall(PetscNew(&chol));

297:   chol->Wrap = MatWrapCholmod_SPQR_seqaij;
298:   B->data    = chol;

300:   B->ops->getinfo = MatGetInfo_CHOLMOD;
301:   B->ops->view    = MatView_CHOLMOD;
302:   B->ops->destroy = MatDestroy_CHOLMOD;

304:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_SPQR));
305:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR));

307:   B->factortype   = MAT_FACTOR_QR;
308:   B->assembled    = PETSC_TRUE;
309:   B->preallocated = PETSC_TRUE;

311:   PetscCall(PetscFree(B->solvertype));
312:   PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype));
313:   B->canuseordering = PETSC_FALSE;
314:   PetscCall(CholmodStart(B));
315:   chol->common->itype = CHOLMOD_LONG;
316:   *F                  = B;
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }