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: /*
10: If A is a MATNORMALHERMITIAN or MATNORMAL format that it pulls the underlying matrix out
11: and performs the factorization on that matrix
12: */
13: static PetscErrorCode MatWrapCholmod_SPQR_seqaij(Mat A, PetscBool values, cholmod_sparse *C, PetscBool *aijalloc, PetscBool *valloc)
14: {
15: Mat_SeqAIJ *aij;
16: Mat AT, B = NULL;
17: Vec left, right;
18: const PetscScalar *aa, *L = NULL;
19: PetscScalar *ca, scale;
20: const PetscInt *ai, *aj;
21: PetscInt n = A->cmap->n, i, j, k, nz;
22: SuiteSparse_long *ci, *cj; /* SuiteSparse_long is the only choice for SPQR */
23: PetscBool vain = PETSC_FALSE, flg;
25: PetscFunctionBegin;
26: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg));
27: if (flg) {
28: PetscCall(MatNormalHermitianGetMat(A, &AT));
29: } else if (!PetscDefined(USE_COMPLEX)) {
30: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg));
31: if (flg) PetscCall(MatNormalGetMat(A, &AT));
32: }
33: if (flg) {
34: B = A;
35: A = AT;
36: 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));
37: PetscCheck(left == right, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalScale() has been called on the input Mat with L != R");
38: if (values && left) {
39: PetscCall(VecGetArrayRead(left, &L));
40: #if PetscDefined(USE_COMPLEX)
41: for (j = 0; j < n; j++)
42: 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");
43: #endif
44: }
45: }
46: /* cholmod_sparse is compressed sparse column */
47: PetscCall(MatIsSymmetric(A, 0.0, &flg));
48: if (flg) {
49: PetscCall(PetscObjectReference((PetscObject)A));
50: AT = A;
51: } else {
52: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
53: }
54: aij = (Mat_SeqAIJ *)AT->data;
55: ai = aij->j;
56: aj = aij->i;
57: for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
58: PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci));
59: if (values) {
60: vain = PETSC_TRUE;
61: PetscCall(PetscMalloc1(nz, &ca));
62: PetscCall(MatSeqAIJGetArrayRead(AT, &aa));
63: }
64: for (j = 0, k = 0; j < n; j++) {
65: cj[j] = k;
66: for (i = aj[j]; i < aj[j + 1]; i++, k++) {
67: ci[k] = ai[i];
68: if (values) {
69: ca[k] = aa[i];
70: if (L) ca[k] *= L[j];
71: }
72: }
73: }
74: cj[j] = k;
75: *aijalloc = PETSC_TRUE;
76: *valloc = vain;
77: if (values) {
78: PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa));
79: if (L) PetscCall(VecRestoreArrayRead(left, &L));
80: }
82: PetscCall(PetscMemzero(C, sizeof(*C)));
84: C->nrow = (size_t)AT->cmap->n;
85: C->ncol = (size_t)AT->rmap->n;
86: C->nzmax = (size_t)nz;
87: C->p = cj;
88: C->i = ci;
89: C->x = values ? ca : 0;
90: C->stype = 0;
91: C->itype = CHOLMOD_LONG;
92: C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
93: C->dtype = CHOLMOD_DOUBLE;
94: C->sorted = 1;
95: C->packed = 1;
97: PetscCall(MatDestroy(&AT));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
102: {
103: PetscFunctionBegin;
104: *type = MATSOLVERSPQR;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: #define GET_ARRAY_READ 0
109: #define GET_ARRAY_WRITE 1
111: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
112: {
113: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
114: cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;
116: PetscFunctionBegin;
117: if (!chol->normal) {
118: QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
119: PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
120: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
121: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
122: } else {
123: Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
124: PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
125: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
126: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
127: PetscCallExternal(!cholmod_l_free_dense, &Z_handle, chol->common);
128: }
129: *_Y_handle = Y_handle;
130: PetscCallExternal(!cholmod_l_free_dense, &QTB_handle, chol->common);
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X)
135: {
136: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
137: cholmod_dense cholB, *Y_handle = NULL;
138: PetscInt n;
139: PetscScalar *v;
141: PetscFunctionBegin;
142: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
143: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
144: PetscCall(VecGetLocalSize(X, &n));
145: PetscCall(VecGetArrayWrite(X, &v));
146: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
147: PetscCall(VecRestoreArrayWrite(X, &v));
148: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
149: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
150: PetscCall(VecScale(X, chol->scale));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
155: {
156: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
157: cholmod_dense cholB, *Y_handle = NULL;
158: PetscScalar *v;
159: PetscInt lda;
161: PetscFunctionBegin;
162: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
163: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
164: PetscCall(MatDenseGetArrayWrite(X, &v));
165: PetscCall(MatDenseGetLDA(X, &lda));
166: if ((size_t)lda == Y_handle->d) {
167: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
168: } else {
169: 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));
170: }
171: PetscCall(MatDenseRestoreArrayWrite(X, &v));
172: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
173: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
174: PetscCall(MatScale(X, chol->scale));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
179: {
180: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
181: cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;
183: PetscFunctionBegin;
184: RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
185: PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
186: Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
187: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
188: *_Y_handle = Y_handle;
189: PetscCallExternal(!cholmod_l_free_dense, &RTB_handle, chol->common);
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
194: {
195: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
196: cholmod_dense cholB, *Y_handle = NULL;
197: PetscInt n;
198: PetscScalar *v;
200: PetscFunctionBegin;
201: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
202: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
203: PetscCall(VecGetLocalSize(X, &n));
204: PetscCall(VecGetArrayWrite(X, &v));
205: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
206: PetscCall(VecRestoreArrayWrite(X, &v));
207: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
208: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
213: {
214: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
215: cholmod_dense cholB, *Y_handle = NULL;
216: PetscScalar *v;
217: PetscInt lda;
219: PetscFunctionBegin;
220: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
221: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
222: PetscCall(MatDenseGetArrayWrite(X, &v));
223: PetscCall(MatDenseGetLDA(X, &lda));
224: if ((size_t)lda == Y_handle->d) {
225: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
226: } else {
227: 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));
228: }
229: PetscCall(MatDenseRestoreArrayWrite(X, &v));
230: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
231: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
236: {
237: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
238: cholmod_sparse cholA;
239: PetscBool aijalloc, valloc;
240: int err;
242: PetscFunctionBegin;
243: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
244: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
245: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
246: err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);
247: PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status);
249: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
250: if (valloc) PetscCall(PetscFree(cholA.x));
252: F->ops->solve = MatSolve_SPQR;
253: F->ops->matsolve = MatMatSolve_SPQR;
254: if (chol->normal) {
255: PetscScalar scale;
257: 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));
258: F->ops->solvetranspose = MatSolve_SPQR;
259: F->ops->matsolvetranspose = MatMatSolve_SPQR;
260: chol->scale = 1.0 / scale;
261: } else if (A->cmap->n == A->rmap->n) {
262: F->ops->solvetranspose = MatSolveTranspose_SPQR;
263: F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR;
264: chol->scale = 1.0;
265: }
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
270: {
271: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
272: cholmod_sparse cholA;
273: PetscBool aijalloc, valloc;
275: PetscFunctionBegin;
276: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
277: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
278: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
279: if (PetscDefined(USE_DEBUG)) PetscCallExternal(!cholmod_l_check_sparse, &cholA, chol->common);
280: if (chol->spqrfact) PetscCallExternal(!SuiteSparseQR_C_free, &chol->spqrfact, chol->common);
281: chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common);
282: PetscCheck(chol->spqrfact, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status);
284: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
285: if (valloc) PetscCall(PetscFree(cholA.x));
287: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*MC
292: MATSOLVERSPQR
294: A matrix solver type providing direct solvers, QR factorizations, for sequential `MATSEQAIJ` and `MATNORMAL` matrices (that contain a `MATSEQAIJ`).
295: via the external package SPQR.
297: Use `./configure --download-suitesparse` to install PETSc to use SPQR
299: Level: beginner
301: Note:
302: SPQR is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>
304: .seealso: [](ch_matrices), `Mat`, `PCQR`, `PCFactorSetMatSolverType()`, `MatSolverType`
305: M*/
307: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A, MatFactorType ftype, Mat *F)
308: {
309: Mat B;
310: Mat_CHOLMOD *chol;
311: PetscInt m = A->rmap->n, n = A->cmap->n;
312: const char *prefix;
314: PetscFunctionBegin;
315: /* Create the factorization matrix F */
316: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
317: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
318: PetscCall(PetscStrallocpy("spqr", &((PetscObject)B)->type_name));
319: PetscCall(MatGetOptionsPrefix(A, &prefix));
320: PetscCall(MatSetOptionsPrefix(B, prefix));
321: PetscCall(MatSetUp(B));
322: PetscCall(PetscNew(&chol));
324: chol->Wrap = MatWrapCholmod_SPQR_seqaij;
325: B->data = chol;
327: B->ops->getinfo = MatGetInfo_CHOLMOD;
328: B->ops->view = MatView_CHOLMOD;
329: B->ops->destroy = MatDestroy_CHOLMOD;
331: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_SPQR));
332: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR));
334: B->factortype = MAT_FACTOR_QR;
335: B->assembled = PETSC_TRUE;
336: B->preallocated = PETSC_TRUE;
338: PetscCall(PetscFree(B->solvertype));
339: PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype));
340: B->canuseordering = PETSC_FALSE;
341: PetscCall(CholmodStart(B));
342: chol->common->itype = CHOLMOD_LONG;
343: *F = B;
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }