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