Actual source code: svdscalap.c

slepc-3.23.1 2025-05-01
Report Typos and Errors
  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   info,lwork=-1,one=1;
 51:   PetscInt       M,N,m,n,mn;
 52: #if defined(PETSC_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 !defined(PETSC_USE_COMPLEX)
 76:   /* allocate workspace */
 77:   PetscCallBLAS("SCALAPACKgesvd",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:   PetscCheckScaLapackInfo("gesvd",info);
 79:   PetscCall(PetscBLASIntCast((PetscInt)minlwork,&lwork));
 80:   PetscCall(PetscMalloc1(lwork,&work));
 81:   /* call computational routine */
 82:   PetscCallBLAS("SCALAPACKgesvd",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));
 83:   PetscCheckScaLapackInfo("gesvd",info);
 84:   PetscCall(PetscFree(work));
 85: #else
 86:   /* allocate workspace */
 87:   PetscCallBLAS("SCALAPACKgesvd",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));
 88:   PetscCheckScaLapackInfo("gesvd",info);
 89:   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(minlwork),&lwork));
 90:   lrwork = 1+4*PetscMax(a->M,a->N);
 91:   PetscCall(PetscMalloc2(lwork,&work,lrwork,&rwork));
 92:   /* call computational routine */
 93:   PetscCallBLAS("SCALAPACKgesvd",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));
 94:   PetscCheckScaLapackInfo("gesvd",info);
 95:   PetscCall(PetscFree2(work,rwork));
 96: #endif
 97:   PetscCall(PetscFPTrapPop());

 99:   PetscCall(MatHermitianTranspose(QT,MAT_INITIAL_MATRIX,&Q));
100:   PetscCall(MatDestroy(&QT));
101:   PetscCall(BVGetMat(svd->U,&U));
102:   PetscCall(BVGetMat(svd->V,&V));
103:   if (M>=N) {
104:     PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&U));
105:     PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&V));
106:   } else {
107:     PetscCall(MatConvert(Q,MATDENSE,MAT_REUSE_MATRIX,&U));
108:     PetscCall(MatConvert(Z,MATDENSE,MAT_REUSE_MATRIX,&V));
109:   }
110:   PetscCall(BVRestoreMat(svd->U,&U));
111:   PetscCall(BVRestoreMat(svd->V,&V));
112:   PetscCall(MatDestroy(&Z));
113:   PetscCall(MatDestroy(&Q));

115:   svd->nconv  = svd->ncv;
116:   svd->its    = 1;
117:   svd->reason = SVD_CONVERGED_TOL;
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: static PetscErrorCode SVDDestroy_ScaLAPACK(SVD svd)
122: {
123:   PetscFunctionBegin;
124:   PetscCall(PetscFree(svd->data));
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: static PetscErrorCode SVDReset_ScaLAPACK(SVD svd)
129: {
130:   SVD_ScaLAPACK  *ctx = (SVD_ScaLAPACK*)svd->data;

132:   PetscFunctionBegin;
133:   PetscCall(MatDestroy(&ctx->As));
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

137: SLEPC_EXTERN PetscErrorCode SVDCreate_ScaLAPACK(SVD svd)
138: {
139:   SVD_ScaLAPACK  *ctx;

141:   PetscFunctionBegin;
142:   PetscCall(PetscNew(&ctx));
143:   svd->data = (void*)ctx;

145:   svd->ops->solve          = SVDSolve_ScaLAPACK;
146:   svd->ops->setup          = SVDSetUp_ScaLAPACK;
147:   svd->ops->destroy        = SVDDestroy_ScaLAPACK;
148:   svd->ops->reset          = SVDReset_ScaLAPACK;
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }