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;
13: const PetscScalar *aa;
14: PetscScalar *ca;
15: const PetscInt *ai, *aj;
16: PetscInt n = A->cmap->n, i, j, k, nz;
17: SuiteSparse_long *ci, *cj; /* SuiteSparse_long is the only choice for SPQR */
18: PetscBool vain = PETSC_FALSE, flg;
20: PetscFunctionBegin;
21: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &flg));
22: if (flg) {
23: PetscCall(MatNormalHermitianGetMat(A, &A));
24: } else if (!PetscDefined(USE_COMPLEX)) {
25: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &flg));
26: if (flg) PetscCall(MatNormalGetMat(A, &A));
27: }
28: /* cholmod_sparse is compressed sparse column */
29: PetscCall(MatIsSymmetric(A, 0.0, &flg));
30: if (flg) {
31: PetscCall(PetscObjectReference((PetscObject)A));
32: AT = A;
33: } else {
34: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
35: }
36: aij = (Mat_SeqAIJ *)AT->data;
37: ai = aij->j;
38: aj = aij->i;
39: for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
40: PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci));
41: if (values) {
42: vain = PETSC_TRUE;
43: PetscCall(PetscMalloc1(nz, &ca));
44: PetscCall(MatSeqAIJGetArrayRead(AT, &aa));
45: }
46: for (j = 0, k = 0; j < n; j++) {
47: cj[j] = k;
48: for (i = aj[j]; i < aj[j + 1]; i++, k++) {
49: ci[k] = ai[i];
50: if (values) ca[k] = aa[i];
51: }
52: }
53: cj[j] = k;
54: *aijalloc = PETSC_TRUE;
55: *valloc = vain;
56: if (values) PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa));
58: PetscCall(PetscMemzero(C, sizeof(*C)));
60: C->nrow = (size_t)AT->cmap->n;
61: C->ncol = (size_t)AT->rmap->n;
62: C->nzmax = (size_t)nz;
63: C->p = cj;
64: C->i = ci;
65: C->x = values ? ca : 0;
66: C->stype = 0;
67: C->itype = CHOLMOD_LONG;
68: C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
69: C->dtype = CHOLMOD_DOUBLE;
70: C->sorted = 1;
71: C->packed = 1;
73: PetscCall(MatDestroy(&AT));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
78: {
79: PetscFunctionBegin;
80: *type = MATSOLVERSPQR;
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: #define GET_ARRAY_READ 0
85: #define GET_ARRAY_WRITE 1
87: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
88: {
89: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
90: cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;
92: PetscFunctionBegin;
93: if (!chol->normal) {
94: QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
95: PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
96: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
97: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
98: } else {
99: Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
100: PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
101: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
102: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
103: PetscCallExternal(!cholmod_l_free_dense, &Z_handle, chol->common);
104: }
105: *_Y_handle = Y_handle;
106: PetscCallExternal(!cholmod_l_free_dense, &QTB_handle, chol->common);
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X)
111: {
112: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
113: cholmod_dense cholB, *Y_handle = NULL;
114: PetscInt n;
115: PetscScalar *v;
117: PetscFunctionBegin;
118: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
119: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
120: PetscCall(VecGetLocalSize(X, &n));
121: PetscCall(VecGetArrayWrite(X, &v));
122: PetscCall(PetscArraycpy(v, (PetscScalar *)(Y_handle->x), n));
123: PetscCall(VecRestoreArrayWrite(X, &v));
124: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
125: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
130: {
131: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
132: cholmod_dense cholB, *Y_handle = NULL;
133: PetscScalar *v;
134: PetscInt lda;
136: PetscFunctionBegin;
137: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
138: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
139: PetscCall(MatDenseGetArrayWrite(X, &v));
140: PetscCall(MatDenseGetLDA(X, &lda));
141: if ((size_t)lda == Y_handle->d) {
142: PetscCall(PetscArraycpy(v, (PetscScalar *)(Y_handle->x), lda * Y_handle->ncol));
143: } else {
144: 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));
145: }
146: PetscCall(MatDenseRestoreArrayWrite(X, &v));
147: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
148: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
153: {
154: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
155: cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;
157: PetscFunctionBegin;
158: RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
159: PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
160: Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
161: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
162: *_Y_handle = Y_handle;
163: PetscCallExternal(!cholmod_l_free_dense, &RTB_handle, chol->common);
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
168: {
169: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
170: cholmod_dense cholB, *Y_handle = NULL;
171: PetscInt n;
172: PetscScalar *v;
174: PetscFunctionBegin;
175: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
176: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
177: PetscCall(VecGetLocalSize(X, &n));
178: PetscCall(VecGetArrayWrite(X, &v));
179: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
180: PetscCall(VecRestoreArrayWrite(X, &v));
181: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
182: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
187: {
188: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
189: cholmod_dense cholB, *Y_handle = NULL;
190: PetscScalar *v;
191: PetscInt lda;
193: PetscFunctionBegin;
194: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
195: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
196: PetscCall(MatDenseGetArrayWrite(X, &v));
197: PetscCall(MatDenseGetLDA(X, &lda));
198: if ((size_t)lda == Y_handle->d) {
199: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
200: } else {
201: 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));
202: }
203: PetscCall(MatDenseRestoreArrayWrite(X, &v));
204: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
205: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
210: {
211: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
212: cholmod_sparse cholA;
213: PetscBool aijalloc, valloc;
214: int err;
216: PetscFunctionBegin;
217: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
218: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
219: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
220: err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);
221: PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status);
223: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
224: if (valloc) PetscCall(PetscFree(cholA.x));
226: F->ops->solve = MatSolve_SPQR;
227: F->ops->matsolve = MatMatSolve_SPQR;
228: if (chol->normal) {
229: F->ops->solvetranspose = MatSolve_SPQR;
230: F->ops->matsolvetranspose = MatMatSolve_SPQR;
231: } else if (A->cmap->n == A->rmap->n) {
232: F->ops->solvetranspose = MatSolveTranspose_SPQR;
233: F->ops->matsolvetranspose = MatMatSolveTranspose_SPQR;
234: }
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: PETSC_INTERN PetscErrorCode MatQRFactorSymbolic_SPQR(Mat F, Mat A, IS perm, const MatFactorInfo *info)
239: {
240: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
241: cholmod_sparse cholA;
242: PetscBool aijalloc, valloc;
244: PetscFunctionBegin;
245: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
246: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
247: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
248: if (PetscDefined(USE_DEBUG)) PetscCallExternal(!cholmod_l_check_sparse, &cholA, chol->common);
249: if (chol->spqrfact) PetscCallExternal(!SuiteSparseQR_C_free, &chol->spqrfact, chol->common);
250: chol->spqrfact = SuiteSparseQR_C_symbolic(SPQR_ORDERING_DEFAULT, 1, &cholA, chol->common);
251: PetscCheck(chol->spqrfact, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status);
253: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
254: if (valloc) PetscCall(PetscFree(cholA.x));
256: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", MatQRFactorNumeric_SPQR));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*MC
261: MATSOLVERSPQR
263: A matrix type providing direct solvers, QR factorizations, for sequential matrices
264: via the external package SPQR.
266: Use `./configure --download-suitesparse` to install PETSc to use SPQR
268: Consult SPQR documentation for more information about the common parameters
269: which correspond to the options database keys below.
271: Level: beginner
273: Note:
274: SPQR is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>
276: .seealso: [](ch_matrices), `Mat`, `PCQR`, `PCFactorSetMatSolverType()`, `MatSolverType`
277: M*/
279: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat A, MatFactorType ftype, Mat *F)
280: {
281: Mat B;
282: Mat_CHOLMOD *chol;
283: PetscInt m = A->rmap->n, n = A->cmap->n;
284: const char *prefix;
286: PetscFunctionBegin;
287: /* Create the factorization matrix F */
288: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
289: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
290: PetscCall(PetscStrallocpy("spqr", &((PetscObject)B)->type_name));
291: PetscCall(MatGetOptionsPrefix(A, &prefix));
292: PetscCall(MatSetOptionsPrefix(B, prefix));
293: PetscCall(MatSetUp(B));
294: PetscCall(PetscNew(&chol));
296: chol->Wrap = MatWrapCholmod_SPQR_seqaij;
297: B->data = chol;
299: B->ops->getinfo = MatGetInfo_CHOLMOD;
300: B->ops->view = MatView_CHOLMOD;
301: B->ops->destroy = MatDestroy_CHOLMOD;
303: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_SPQR));
304: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SPQR));
306: B->factortype = MAT_FACTOR_QR;
307: B->assembled = PETSC_TRUE;
308: B->preallocated = PETSC_TRUE;
310: PetscCall(PetscFree(B->solvertype));
311: PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype));
312: B->canuseordering = PETSC_FALSE;
313: PetscCall(CholmodStart(B));
314: chol->common->itype = CHOLMOD_LONG;
315: *F = B;
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }