Actual source code: aijspqr.c

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

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

  9: static PetscErrorCode MatWrapCholmod_SPQR_seqaij(Mat A, PetscBool values, cholmod_sparse *C, PetscBool *aijalloc, PetscBool *valloc)
 10: {
 11:   Mat_SeqAIJ        *aij;
 12:   Mat                AT, B = NULL;
 13:   Vec                left, right;
 14:   const PetscScalar *aa, *L = NULL;
 15:   PetscScalar       *ca, scale;
 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, &AT));
 25:   } else if (!PetscDefined(USE_COMPLEX)) {
 26:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg));
 27:     if (flg) PetscCall(MatNormalGetMat(A, &AT));
 28:   }
 29:   if (flg) {
 30:     B = A;
 31:     A = AT;
 32:     PetscCall(MatShellGetScalingShifts(B, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, &left, &right, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
 33:     PetscCheck(left == right, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalScale() has been called on the input Mat with L != R");
 34:     if (values && left) {
 35:       PetscCall(VecGetArrayRead(left, &L));
 36: #if PetscDefined(USE_COMPLEX)
 37:       for (j = 0; j < n; j++)
 38:         PetscCheck(PetscAbsReal(PetscImaginaryPart(L[j])) < PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalScale() has been called on the input Mat with a complex Vec");
 39: #endif
 40:     }
 41:   }
 42:   /* cholmod_sparse is compressed sparse column */
 43:   PetscCall(MatIsSymmetric(A, 0.0, &flg));
 44:   if (flg) {
 45:     PetscCall(PetscObjectReference((PetscObject)A));
 46:     AT = A;
 47:   } else {
 48:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
 49:   }
 50:   aij = (Mat_SeqAIJ *)AT->data;
 51:   ai  = aij->j;
 52:   aj  = aij->i;
 53:   for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
 54:   PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci));
 55:   if (values) {
 56:     vain = PETSC_TRUE;
 57:     PetscCall(PetscMalloc1(nz, &ca));
 58:     PetscCall(MatSeqAIJGetArrayRead(AT, &aa));
 59:   }
 60:   for (j = 0, k = 0; j < n; j++) {
 61:     cj[j] = k;
 62:     for (i = aj[j]; i < aj[j + 1]; i++, k++) {
 63:       ci[k] = ai[i];
 64:       if (values) {
 65:         ca[k] = aa[i];
 66:         if (L) ca[k] *= L[j];
 67:       }
 68:     }
 69:   }
 70:   cj[j]     = k;
 71:   *aijalloc = PETSC_TRUE;
 72:   *valloc   = vain;
 73:   if (values) {
 74:     PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa));
 75:     if (L) PetscCall(VecRestoreArrayRead(left, &L));
 76:   }

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

 80:   C->nrow   = (size_t)AT->cmap->n;
 81:   C->ncol   = (size_t)AT->rmap->n;
 82:   C->nzmax  = (size_t)nz;
 83:   C->p      = cj;
 84:   C->i      = ci;
 85:   C->x      = values ? ca : 0;
 86:   C->stype  = 0;
 87:   C->itype  = CHOLMOD_LONG;
 88:   C->xtype  = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
 89:   C->dtype  = CHOLMOD_DOUBLE;
 90:   C->sorted = 1;
 91:   C->packed = 1;

 93:   PetscCall(MatDestroy(&AT));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
 98: {
 99:   PetscFunctionBegin;
100:   *type = MATSOLVERSPQR;
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: #define GET_ARRAY_READ  0
105: #define GET_ARRAY_WRITE 1

107: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
108: {
109:   Mat_CHOLMOD   *chol     = (Mat_CHOLMOD *)F->data;
110:   cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;

112:   PetscFunctionBegin;
113:   if (!chol->normal) {
114:     QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
115:     PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
116:     Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
117:     PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
118:   } else {
119:     Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
120:     PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
121:     Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
122:     PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
123:     PetscCallExternal(!cholmod_l_free_dense, &Z_handle, chol->common);
124:   }
125:   *_Y_handle = Y_handle;
126:   PetscCallExternal(!cholmod_l_free_dense, &QTB_handle, chol->common);
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

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

137:   PetscFunctionBegin;
138:   PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
139:   PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
140:   PetscCall(VecGetLocalSize(X, &n));
141:   PetscCall(VecGetArrayWrite(X, &v));
142:   PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
143:   PetscCall(VecRestoreArrayWrite(X, &v));
144:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
145:   PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
146:   PetscCall(VecScale(X, chol->scale));
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
151: {
152:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
153:   cholmod_dense cholB, *Y_handle = NULL;
154:   PetscScalar  *v;
155:   PetscInt      lda;

157:   PetscFunctionBegin;
158:   PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
159:   PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
160:   PetscCall(MatDenseGetArrayWrite(X, &v));
161:   PetscCall(MatDenseGetLDA(X, &lda));
162:   if ((size_t)lda == Y_handle->d) {
163:     PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
164:   } else {
165:     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));
166:   }
167:   PetscCall(MatDenseRestoreArrayWrite(X, &v));
168:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
169:   PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
170:   PetscCall(MatScale(X, chol->scale));
171:   PetscFunctionReturn(PETSC_SUCCESS);
172: }

174: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
175: {
176:   Mat_CHOLMOD   *chol     = (Mat_CHOLMOD *)F->data;
177:   cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;

179:   PetscFunctionBegin;
180:   RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
181:   PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
182:   Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
183:   PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
184:   *_Y_handle = Y_handle;
185:   PetscCallExternal(!cholmod_l_free_dense, &RTB_handle, chol->common);
186:   PetscFunctionReturn(PETSC_SUCCESS);
187: }

189: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
190: {
191:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
192:   cholmod_dense cholB, *Y_handle = NULL;
193:   PetscInt      n;
194:   PetscScalar  *v;

196:   PetscFunctionBegin;
197:   PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
198:   PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
199:   PetscCall(VecGetLocalSize(X, &n));
200:   PetscCall(VecGetArrayWrite(X, &v));
201:   PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
202:   PetscCall(VecRestoreArrayWrite(X, &v));
203:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
204:   PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
209: {
210:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
211:   cholmod_dense cholB, *Y_handle = NULL;
212:   PetscScalar  *v;
213:   PetscInt      lda;

215:   PetscFunctionBegin;
216:   PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
217:   PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
218:   PetscCall(MatDenseGetArrayWrite(X, &v));
219:   PetscCall(MatDenseGetLDA(X, &lda));
220:   if ((size_t)lda == Y_handle->d) {
221:     PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
222:   } else {
223:     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));
224:   }
225:   PetscCall(MatDenseRestoreArrayWrite(X, &v));
226:   PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
227:   PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
232: {
233:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
234:   cholmod_sparse cholA;
235:   PetscBool      aijalloc, valloc;
236:   int            err;

238:   PetscFunctionBegin;
239:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
240:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
241:   PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
242:   err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);
243:   PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status);

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

248:   F->ops->solve    = MatSolve_SPQR;
249:   F->ops->matsolve = MatMatSolve_SPQR;
250:   if (chol->normal) {
251:     PetscScalar scale;

253:     PetscCall(MatShellGetScalingShifts(A, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, NULL, NULL, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
254:     F->ops->solvetranspose    = MatSolve_SPQR;
255:     F->ops->matsolvetranspose = MatMatSolve_SPQR;
256:     chol->scale               = 1.0 / scale;
257:   } else if (A->cmap->n == A->rmap->n) {
258:     F->ops->solvetranspose    = MatSolveTranspose_SPQR;
259:     F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR;
260:     chol->scale               = 1.0;
261:   }
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

265: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
266: {
267:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
268:   cholmod_sparse cholA;
269:   PetscBool      aijalloc, valloc;

271:   PetscFunctionBegin;
272:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
273:   if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
274:   PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
275:   if (PetscDefined(USE_DEBUG)) PetscCallExternal(!cholmod_l_check_sparse, &cholA, chol->common);
276:   if (chol->spqrfact) PetscCallExternal(!SuiteSparseQR_C_free, &chol->spqrfact, chol->common);
277:   chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common);
278:   PetscCheck(chol->spqrfact, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status);

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

283:   PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR));
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: /*MC
288:   MATSOLVERSPQR

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

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

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

298:    Level: beginner

300:    Note:
301:    SPQR is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>

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

306: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A, MatFactorType ftype, Mat *F)
307: {
308:   Mat          B;
309:   Mat_CHOLMOD *chol;
310:   PetscInt     m = A->rmap->n, n = A->cmap->n;
311:   const char  *prefix;

313:   PetscFunctionBegin;
314:   /* Create the factorization matrix F */
315:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
316:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
317:   PetscCall(PetscStrallocpy("spqr", &((PetscObject)B)->type_name));
318:   PetscCall(MatGetOptionsPrefix(A, &prefix));
319:   PetscCall(MatSetOptionsPrefix(B, prefix));
320:   PetscCall(MatSetUp(B));
321:   PetscCall(PetscNew(&chol));

323:   chol->Wrap = MatWrapCholmod_SPQR_seqaij;
324:   B->data    = chol;

326:   B->ops->getinfo = MatGetInfo_CHOLMOD;
327:   B->ops->view    = MatView_CHOLMOD;
328:   B->ops->destroy = MatDestroy_CHOLMOD;

330:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_SPQR));
331:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR));

333:   B->factortype   = MAT_FACTOR_QR;
334:   B->assembled    = PETSC_TRUE;
335:   B->preallocated = PETSC_TRUE;

337:   PetscCall(PetscFree(B->solvertype));
338:   PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype));
339:   B->canuseordering = PETSC_FALSE;
340:   PetscCall(CholmodStart(B));
341:   chol->common->itype = CHOLMOD_LONG;
342:   *F                  = B;
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }