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, B = NULL;
13: Vec left, right;
14: const PetscScalar *aa, *L = NULL;
15: PetscScalar *ca, scale;
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: 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));
33: PetscCheck(left == right, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Cannot call SuiteSparseQR if MatDiagonalScale() has been called on the input Mat with L != R");
34: if (values && left) {
35: PetscCall(VecGetArrayRead(left, &L));
36: #if PetscDefined(USE_COMPLEX)
37: for (j = 0; j < n; j++)
38: 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");
39: #endif
40: }
41: }
42: /* cholmod_sparse is compressed sparse column */
43: PetscCall(MatIsSymmetric(A, 0.0, &flg));
44: if (flg) {
45: PetscCall(PetscObjectReference((PetscObject)A));
46: AT = A;
47: } else {
48: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
49: }
50: aij = (Mat_SeqAIJ *)AT->data;
51: ai = aij->j;
52: aj = aij->i;
53: for (j = 0, nz = 0; j < n; j++) nz += aj[j + 1] - aj[j];
54: PetscCall(PetscMalloc2(n + 1, &cj, nz, &ci));
55: if (values) {
56: vain = PETSC_TRUE;
57: PetscCall(PetscMalloc1(nz, &ca));
58: PetscCall(MatSeqAIJGetArrayRead(AT, &aa));
59: }
60: for (j = 0, k = 0; j < n; j++) {
61: cj[j] = k;
62: for (i = aj[j]; i < aj[j + 1]; i++, k++) {
63: ci[k] = ai[i];
64: if (values) {
65: ca[k] = aa[i];
66: if (L) ca[k] *= L[j];
67: }
68: }
69: }
70: cj[j] = k;
71: *aijalloc = PETSC_TRUE;
72: *valloc = vain;
73: if (values) {
74: PetscCall(MatSeqAIJRestoreArrayRead(AT, &aa));
75: if (L) PetscCall(VecRestoreArrayRead(left, &L));
76: }
78: PetscCall(PetscMemzero(C, sizeof(*C)));
80: C->nrow = (size_t)AT->cmap->n;
81: C->ncol = (size_t)AT->rmap->n;
82: C->nzmax = (size_t)nz;
83: C->p = cj;
84: C->i = ci;
85: C->x = values ? ca : 0;
86: C->stype = 0;
87: C->itype = CHOLMOD_LONG;
88: C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
89: C->dtype = CHOLMOD_DOUBLE;
90: C->sorted = 1;
91: C->packed = 1;
93: PetscCall(MatDestroy(&AT));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode MatFactorGetSolverType_seqaij_SPQR(Mat A, MatSolverType *type)
98: {
99: PetscFunctionBegin;
100: *type = MATSOLVERSPQR;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: #define GET_ARRAY_READ 0
105: #define GET_ARRAY_WRITE 1
107: static PetscErrorCode MatSolve_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
108: {
109: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
110: cholmod_dense *Y_handle = NULL, *QTB_handle = NULL, *Z_handle = NULL;
112: PetscFunctionBegin;
113: if (!chol->normal) {
114: QTB_handle = SuiteSparseQR_C_qmult(SPQR_QTX, chol->spqrfact, cholB, chol->common);
115: PetscCheck(QTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
116: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, QTB_handle, chol->common);
117: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
118: } else {
119: Z_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
120: PetscCheck(Z_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
121: Y_handle = SuiteSparseQR_C_solve(SPQR_RETX_EQUALS_B, chol->spqrfact, Z_handle, chol->common);
122: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
123: PetscCallExternal(!cholmod_l_free_dense, &Z_handle, chol->common);
124: }
125: *_Y_handle = Y_handle;
126: PetscCallExternal(!cholmod_l_free_dense, &QTB_handle, chol->common);
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode MatSolve_SPQR(Mat F, Vec B, Vec X)
131: {
132: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
133: cholmod_dense cholB, *Y_handle = NULL;
134: PetscInt n;
135: PetscScalar *v;
137: PetscFunctionBegin;
138: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
139: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
140: PetscCall(VecGetLocalSize(X, &n));
141: PetscCall(VecGetArrayWrite(X, &v));
142: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
143: PetscCall(VecRestoreArrayWrite(X, &v));
144: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
145: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
146: PetscCall(VecScale(X, chol->scale));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode MatMatSolve_SPQR(Mat F, Mat B, Mat X)
151: {
152: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
153: cholmod_dense cholB, *Y_handle = NULL;
154: PetscScalar *v;
155: PetscInt lda;
157: PetscFunctionBegin;
158: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
159: PetscCall(MatSolve_SPQR_Internal(F, &cholB, &Y_handle));
160: PetscCall(MatDenseGetArrayWrite(X, &v));
161: PetscCall(MatDenseGetLDA(X, &lda));
162: if ((size_t)lda == Y_handle->d) {
163: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
164: } else {
165: 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));
166: }
167: PetscCall(MatDenseRestoreArrayWrite(X, &v));
168: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
169: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
170: PetscCall(MatScale(X, chol->scale));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode MatSolveTranspose_SPQR_Internal(Mat F, cholmod_dense *cholB, cholmod_dense **_Y_handle)
175: {
176: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
177: cholmod_dense *Y_handle = NULL, *RTB_handle = NULL;
179: PetscFunctionBegin;
180: RTB_handle = SuiteSparseQR_C_solve(SPQR_RTX_EQUALS_ETB, chol->spqrfact, cholB, chol->common);
181: PetscCheck(RTB_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_solve failed");
182: Y_handle = SuiteSparseQR_C_qmult(SPQR_QX, chol->spqrfact, RTB_handle, chol->common);
183: PetscCheck(Y_handle, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SuiteSparseQR_C_qmult failed");
184: *_Y_handle = Y_handle;
185: PetscCallExternal(!cholmod_l_free_dense, &RTB_handle, chol->common);
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode MatSolveTranspose_SPQR(Mat F, Vec B, Vec X)
190: {
191: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
192: cholmod_dense cholB, *Y_handle = NULL;
193: PetscInt n;
194: PetscScalar *v;
196: PetscFunctionBegin;
197: PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
198: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
199: PetscCall(VecGetLocalSize(X, &n));
200: PetscCall(VecGetArrayWrite(X, &v));
201: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, n));
202: PetscCall(VecRestoreArrayWrite(X, &v));
203: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
204: PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode MatMatSolveTranspose_SPQR(Mat F, Mat B, Mat X)
209: {
210: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
211: cholmod_dense cholB, *Y_handle = NULL;
212: PetscScalar *v;
213: PetscInt lda;
215: PetscFunctionBegin;
216: PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
217: PetscCall(MatSolveTranspose_SPQR_Internal(F, &cholB, &Y_handle));
218: PetscCall(MatDenseGetArrayWrite(X, &v));
219: PetscCall(MatDenseGetLDA(X, &lda));
220: if ((size_t)lda == Y_handle->d) {
221: PetscCall(PetscArraycpy(v, (PetscScalar *)Y_handle->x, lda * Y_handle->ncol));
222: } else {
223: 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));
224: }
225: PetscCall(MatDenseRestoreArrayWrite(X, &v));
226: PetscCallExternal(!cholmod_l_free_dense, &Y_handle, chol->common);
227: PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode MatQRFactorNumeric_SPQR(Mat F, Mat A, const MatFactorInfo *info)
232: {
233: Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;
234: cholmod_sparse cholA;
235: PetscBool aijalloc, valloc;
236: int err;
238: PetscFunctionBegin;
239: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMALHERMITIAN, &chol->normal));
240: if (!chol->normal && !PetscDefined(USE_COMPLEX)) PetscCall(PetscObjectTypeCompare((PetscObject)A, MATNORMAL, &chol->normal));
241: PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
242: err = !SuiteSparseQR_C_numeric(PETSC_SMALL, &cholA, chol->spqrfact, chol->common);
243: PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "SPQR factorization failed with status %d", chol->common->status);
245: if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
246: if (valloc) PetscCall(PetscFree(cholA.x));
248: F->ops->solve = MatSolve_SPQR;
249: F->ops->matsolve = MatMatSolve_SPQR;
250: if (chol->normal) {
251: PetscScalar scale;
253: 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));
254: F->ops->solvetranspose = MatSolve_SPQR;
255: F->ops->matsolvetranspose = MatMatSolve_SPQR;
256: chol->scale = 1.0 / scale;
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: }