Actual source code: svdscalap.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: This file implements a wrapper to the ScaLAPACK SVD solver
12: */
14: #include <slepc/private/svdimpl.h>
15: #include <slepc/private/slepcscalapack.h>
17: typedef struct {
18: Mat As; /* converted matrix */
19: } SVD_ScaLAPACK;
21: static PetscErrorCode SVDSetUp_ScaLAPACK(SVD svd)
22: {
23: SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
24: PetscInt M,N;
26: PetscFunctionBegin;
27: SVDCheckStandard(svd);
28: SVDCheckDefinite(svd);
29: if (svd->nsv==0) svd->nsv = 1;
30: PetscCall(MatGetSize(svd->A,&M,&N));
31: svd->ncv = N;
32: if (svd->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
33: if (svd->max_it==PETSC_DETERMINE) svd->max_it = 1;
34: svd->leftbasis = PETSC_TRUE;
35: SVDCheckIgnored(svd,SVD_FEATURE_STOPPING);
36: PetscCall(SVDAllocateSolution(svd,0));
38: /* convert matrix */
39: PetscCall(MatDestroy(&ctx->As));
40: PetscCall(MatConvert(svd->OP,MATSCALAPACK,MAT_INITIAL_MATRIX,&ctx->As));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode SVDSolve_ScaLAPACK(SVD svd)
45: {
46: SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
47: Mat A = ctx->As,Z,Q,QT,U,V;
48: Mat_ScaLAPACK *a = (Mat_ScaLAPACK*)A->data,*q,*z;
49: PetscScalar *work,minlwork;
50: PetscBLASInt lwork=-1,one=1;
51: PetscInt M,N,m,n,mn;
52: #if PetscDefined(USE_COMPLEX)
53: PetscBLASInt lrwork;
54: PetscReal *rwork,dummy;
55: #endif
57: PetscFunctionBegin;
58: PetscCall(MatGetSize(A,&M,&N));
59: PetscCall(MatGetLocalSize(A,&m,&n));
60: mn = (M>=N)? n: m;
61: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&Z));
62: PetscCall(MatSetSizes(Z,m,mn,PETSC_DECIDE,PETSC_DECIDE));
63: PetscCall(MatSetType(Z,MATSCALAPACK));
64: PetscCall(MatAssemblyBegin(Z,MAT_FINAL_ASSEMBLY));
65: PetscCall(MatAssemblyEnd(Z,MAT_FINAL_ASSEMBLY));
66: z = (Mat_ScaLAPACK*)Z->data;
67: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&QT));
68: PetscCall(MatSetSizes(QT,mn,n,PETSC_DECIDE,PETSC_DECIDE));
69: PetscCall(MatSetType(QT,MATSCALAPACK));
70: PetscCall(MatAssemblyBegin(QT,MAT_FINAL_ASSEMBLY));
71: PetscCall(MatAssemblyEnd(QT,MAT_FINAL_ASSEMBLY));
72: q = (Mat_ScaLAPACK*)QT->data;
74: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
75: #if !PetscDefined(USE_COMPLEX)
76: /* allocate workspace */
77: PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,&minlwork,&lwork,&info));
78: PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
79: PetscCall(PetscMalloc1(lwork,&work));
80: /* call computational routine */
81: PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,work,&lwork,&info));
82: PetscCall(PetscFree(work));
83: #else
84: /* allocate workspace */
85: PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,&minlwork,&lwork,&dummy,&info));
86: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork),&lwork));
87: lrwork = 1+4*PetscMax(a->M,a->N);
88: PetscCall(PetscMalloc2(lwork,&work,lrwork,&rwork));
89: /* call computational routine */
90: PetscCallScaLAPACKInfo("gesvd",SCALAPACKgesvd_("V","V",&a->M,&a->N,a->loc,&one,&one,a->desc,svd->sigma,z->loc,&one,&one,z->desc,q->loc,&one,&one,q->desc,work,&lwork,rwork,&info));
91: PetscCall(PetscFree2(work,rwork));
92: #endif
93: PetscCall(PetscFPTrapPop());
95: PetscCall(MatHermitianTranspose(QT,MAT_INITIAL_MATRIX,&Q));
96: PetscCall(MatDestroy(&QT));
97: PetscCall(BVGetMat(svd->U,&U));
98: PetscCall(BVGetMat(svd->V,&V));
99: if (M>=N) {
100: PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
101: PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
102: } else {
103: PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&U));
104: PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&V));
105: }
106: PetscCall(BVRestoreMat(svd->U,&U));
107: PetscCall(BVRestoreMat(svd->V,&V));
108: PetscCall(MatDestroy(&Z));
109: PetscCall(MatDestroy(&Q));
111: svd->nconv = svd->ncv;
112: svd->its = 1;
113: svd->reason = SVD_CONVERGED_TOL;
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode SVDDestroy_ScaLAPACK(SVD svd)
118: {
119: PetscFunctionBegin;
120: PetscCall(PetscFree(svd->data));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode SVDReset_ScaLAPACK(SVD svd)
125: {
126: SVD_ScaLAPACK *ctx = (SVD_ScaLAPACK*)svd->data;
128: PetscFunctionBegin;
129: PetscCall(MatDestroy(&ctx->As));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*MC
134: SVDSCALAPACK - SVDSCALAPACK = "scalapack" - A wrapper to the ScaLAPACK
135: singular value solver {cite:p}`Bla97`.
137: Notes:
138: Only available for standard SVD problems, using subroutine `pdgesvd` and
139: the analogs for other precisions.
141: This is a direct singular value solver, that is, the full decomposition
142: is computed. The computation involves redistributing the matrices from PETSc
143: storage to ScaLAPACK distribution, and vice versa (this is done automatically
144: by SLEPc). Alternatively, the user may create the problem matrices
145: already with type `MATSCALAPACK`.
147: Level: beginner
149: .seealso: [](ch:svd), `SVD`, `SVDType`, `SVDSetType()`
150: M*/
151: SLEPC_EXTERN PetscErrorCode SVDCreate_ScaLAPACK(SVD svd)
152: {
153: SVD_ScaLAPACK *ctx;
155: PetscFunctionBegin;
156: PetscCall(PetscNew(&ctx));
157: svd->data = (void*)ctx;
159: svd->ops->solve = SVDSolve_ScaLAPACK;
160: svd->ops->setup = SVDSetUp_ScaLAPACK;
161: svd->ops->destroy = SVDDestroy_ScaLAPACK;
162: svd->ops->reset = SVDReset_ScaLAPACK;
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }