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>
4: #include <../src/mat/impls/shell/shell.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, B = NULL;
14: const PetscScalar *aa, *L = NULL;
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, &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: PetscCheck(!((Mat_Shell *)B->data)->zrows && !((Mat_Shell *)B->data)->zcols, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME
33: PetscCheck(!((Mat_Shell *)B->data)->axpy, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatAXPY() has been called on the input Mat"); // TODO FIXME
34: PetscCheck(((Mat_Shell *)B->data)->left == ((Mat_Shell *)B->data)->right, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalScale() has been called on the input Mat with L != R"); // TODO FIXME
35: if (values && ((Mat_Shell *)B->data)->left) {
36: PetscCall(VecGetArrayRead(((Mat_Shell *)B->data)->left, &L));
37: #if PetscDefined(USE_COMPLEX)
38: for (j = 0; j < n; j++)
39: 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");
40: #endif
41: }
42: PetscCheck(!((Mat_Shell *)B->data)->dshift, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME
43: PetscCheck(PetscAbsScalar(((Mat_Shell *)B->data)->vshift) < PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatShift() has been called on the input Mat"); // TODO FIXME
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(((Mat_Shell *)B->data)->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: F->ops->solvetranspose = MatSolve_SPQR;
255: F->ops->matsolvetranspose = MatMatSolve_SPQR;
256: chol->scale = 1.0 / ((Mat_Shell *)A->data)->vscale;
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: }