Actual source code: aijspqr.c
1: #include <../src/mat/impls/aij/seq/aij.h>
2: #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h>
4: EXTERN_C_BEGIN
5: #include <SuiteSparseQR_C.h>
6: EXTERN_C_END
8: /*
9: If A is a MATNORMALHERMITIAN or MATNORMAL format that it pulls the underlying matrix out
10: and performs the factorization on that matrix
11: */
12: static PetscErrorCode MatWrapCholmod_SPQR_seqaij(Mat A, PetscBool values, cholmod_sparse *C, PetscBool *aijalloc, PetscBool *valloc)
13: {
14: Mat_SeqAIJ *aij;
15: Mat AT, B = NULL;
16: Vec left, right;
17: const PetscScalar *aa, *L = NULL;
18: PetscScalar *ca, scale;
19: const PetscInt *ai, *aj;
20: PetscInt n = A->cmap->n, i, j, k, nz;
21: SuiteSparse_long *ci, *cj; /* SuiteSparse_long is the only choice for SPQR */
22: PetscBool vain = PETSC_FALSE, flg;
24: PetscFunctionBegin;
25: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg));
26: if (flg) {
27: PetscCall(MatNormalHermitianGetMat(A, &AT));
28: } else if (!PetscDefined(USE_COMPLEX)) {
29: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg));
30: if (flg) PetscCall(MatNormalGetMat(A, &AT));
31: }
32: if (flg) {
33: B = A;
34: A = AT;
35: 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));
36: PetscCheck(left == right, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalScale() has been called on the input Mat with L != R");
37: if (values && left) {
38: PetscCall(VecGetArrayRead(left, &L));
39: #if PetscDefined(USE_COMPLEX)
40: for (j = 0; j < n; j++)
41: 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");
42: #endif
43: }
44: }
45: /* cholmod_sparse is compressed sparse column */
46: PetscCall(MatIsSymmetric(A, 0.0, &flg));
47: if (flg) {
48: PetscCall(PetscObjectReference((PetscObject)A));
49: AT = A;
50: } else {
51: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
52: }
53: aij = (Mat_SeqAIJ *)AT->data;
54: ai = aij->j;
55: aj = aij->i;
56: for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
57: PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci));
58: if (values) {
59: vain = PETSC_TRUE;
60: PetscCall(PetscMalloc1(nz, &ca));
61: PetscCall(MatSeqAIJGetArrayRead(AT, &aa));
62: }
63: for (j = 0, k = 0; j < n; j++) {
64: cj[j] = k;
65: for (i = aj[j]; i < aj[j + 1]; i++, k++) {
66: ci[k] = ai[i];
67: if (values) {
68: ca[k] = aa[i];
69: if (L) ca[k] *= L[j];
70: }
71: }
72: }
73: cj[j] = k;
74: *aijalloc = PETSC_TRUE;
75: *valloc = vain;
76: if (values) {
77: PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa));
78: if (L) PetscCall(VecRestoreArrayRead(left, &L));
79: }
81: PetscCall(PetscMemzero(C, sizeof(*C)));
83: C->nrow = (size_t)AT->cmap->n;
84: C->ncol = (size_t)AT->rmap->n;
85: C->nzmax = (size_t)nz;
86: C->p = cj;
87: C->i = ci;
88: C->x = values ? ca : 0;
89: C->stype = 0;
90: C->itype = CHOLMOD_LONG;
91: C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
92: C->dtype = CHOLMOD_DOUBLE;
93: C->sorted = 1;
94: C->packed = 1;
96: PetscCall(MatDestroy(&AT));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
101: {
102: PetscFunctionBegin;
103: *type = MATSOLVERSPQR;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: #define GET_ARRAY_READ 0
108: #define GET_ARRAY_WRITE 1
110: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
111: {
112: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
113: cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;
115: PetscFunctionBegin;
116: if (!chol->normal) {
117: QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
118: PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
119: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
120: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
121: } else {
122: Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
123: PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
124: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
125: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
126: PetscCallExternal(!cholmod_l_free_dense, &Z_handle, chol->common);
127: }
128: *_Y_handle = Y_handle;
129: PetscCallExternal(!cholmod_l_free_dense, &QTB_handle, chol->common);
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X)
134: {
135: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
136: cholmod_dense cholB, *Y_handle = NULL;
137: PetscInt n;
138: PetscScalar *v;
140: PetscFunctionBegin;
141: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
142: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
143: PetscCall(VecGetLocalSize(X, &n));
144: PetscCall(VecGetArrayWrite(X, &v));
145: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
146: PetscCall(VecRestoreArrayWrite(X, &v));
147: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
148: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
149: PetscCall(VecScale(X, chol->scale));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
154: {
155: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
156: cholmod_dense cholB, *Y_handle = NULL;
157: PetscScalar *v;
158: PetscInt lda;
160: PetscFunctionBegin;
161: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
162: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
163: PetscCall(MatDenseGetArrayWrite(X, &v));
164: PetscCall(MatDenseGetLDA(X, &lda));
165: if ((size_t)lda == Y_handle->d) {
166: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
167: } else {
168: 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));
169: }
170: PetscCall(MatDenseRestoreArrayWrite(X, &v));
171: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
172: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
173: PetscCall(MatScale(X, chol->scale));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
178: {
179: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
180: cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;
182: PetscFunctionBegin;
183: RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
184: PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
185: Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
186: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
187: *_Y_handle = Y_handle;
188: PetscCallExternal(!cholmod_l_free_dense, &RTB_handle, chol->common);
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
193: {
194: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
195: cholmod_dense cholB, *Y_handle = NULL;
196: PetscInt n;
197: PetscScalar *v;
199: PetscFunctionBegin;
200: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
201: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
202: PetscCall(VecGetLocalSize(X, &n));
203: PetscCall(VecGetArrayWrite(X, &v));
204: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
205: PetscCall(VecRestoreArrayWrite(X, &v));
206: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
207: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
212: {
213: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
214: cholmod_dense cholB, *Y_handle = NULL;
215: PetscScalar *v;
216: PetscInt lda;
218: PetscFunctionBegin;
219: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
220: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
221: PetscCall(MatDenseGetArrayWrite(X, &v));
222: PetscCall(MatDenseGetLDA(X, &lda));
223: if ((size_t)lda == Y_handle->d) {
224: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
225: } else {
226: 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));
227: }
228: PetscCall(MatDenseRestoreArrayWrite(X, &v));
229: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
230: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
235: {
236: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
237: cholmod_sparse cholA;
238: PetscBool aijalloc, valloc;
239: int err;
241: PetscFunctionBegin;
242: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
243: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
244: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
245: err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);
246: PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status);
248: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
249: if (valloc) PetscCall(PetscFree(cholA.x));
251: F->ops->solve = MatSolve_SPQR;
252: F->ops->matsolve = MatMatSolve_SPQR;
253: if (chol->normal) {
254: PetscScalar scale;
256: 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));
257: F->ops->solvetranspose = MatSolve_SPQR;
258: F->ops->matsolvetranspose = MatMatSolve_SPQR;
259: chol->scale = 1.0 / scale;
260: } else if (A->cmap->n == A->rmap->n) {
261: F->ops->solvetranspose = MatSolveTranspose_SPQR;
262: F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR;
263: chol->scale = 1.0;
264: }
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
269: {
270: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
271: cholmod_sparse cholA;
272: PetscBool aijalloc, valloc;
274: PetscFunctionBegin;
275: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
276: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
277: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
278: if (PetscDefined(USE_DEBUG)) PetscCallExternal(!cholmod_l_check_sparse, &cholA, chol->common);
279: if (chol->spqrfact) PetscCallExternal(!SuiteSparseQR_C_free, &chol->spqrfact, chol->common);
280: chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common);
281: PetscCheck(chol->spqrfact, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status);
283: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
284: if (valloc) PetscCall(PetscFree(cholA.x));
286: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /*MC
291: MATSOLVERSPQR
293: A matrix solver type providing direct solvers, QR factorizations, for sequential `MATSEQAIJ` and `MATNORMAL` matrices (that contain a `MATSEQAIJ`).
294: via the external package SPQR.
296: Use `./configure --download-suitesparse` to install PETSc to use SPQR
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: }