Actual source code: svdlapack.c

slepc-3.23.0 2025-03-29
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 LAPACK SVD subroutines
 12: */

 14: #include <slepc/private/svdimpl.h>
 15: #include <slepcblaslapack.h>

 17: static PetscErrorCode SVDSetUp_LAPACK(SVD svd)
 18: {
 19:   PetscInt       M,N,P=0;

 21:   PetscFunctionBegin;
 22:   if (svd->nsv==0) svd->nsv = 1;
 23:   PetscCall(MatGetSize(svd->A,&M,&N));
 24:   if (!svd->isgeneralized) svd->ncv = N;
 25:   else {
 26:     PetscCall(MatGetSize(svd->OPb,&P,NULL));
 27:     svd->ncv = PetscMin(M,PetscMin(N,P));
 28:   }
 29:   if (svd->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
 30:   SVDCheckIgnored(svd,SVD_FEATURE_STOPPING);
 31:   if (svd->max_it==PETSC_DETERMINE) svd->max_it = 1;
 32:   svd->leftbasis = PETSC_TRUE;
 33:   PetscCall(SVDAllocateSolution(svd,0));
 34:   PetscCall(DSAllocate(svd->ds,PetscMax(N,PetscMax(M,P))));
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

 38: static PetscErrorCode SVDSolve_LAPACK(SVD svd)
 39: {
 40:   PetscInt          M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
 41:   Mat               A,Ar,mat;
 42:   Vec               u,v;
 43:   PetscScalar       *pU,*pV,*pu,*pv,*w;

 45:   PetscFunctionBegin;
 46:   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
 47:   PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
 48:   PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat));
 49:   PetscCall(MatDestroy(&Ar));
 50:   PetscCall(MatGetSize(mat,&M,&N));
 51:   PetscCall(DSSetDimensions(svd->ds,M,0,0));
 52:   PetscCall(DSSVDSetDimensions(svd->ds,N));
 53:   PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
 54:   PetscCall(MatCopy(mat,A,SAME_NONZERO_PATTERN));
 55:   PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
 56:   PetscCall(DSSetState(svd->ds,DS_STATE_RAW));

 58:   n = PetscMin(M,N);
 59:   PetscCall(PetscMalloc1(n,&w));
 60:   PetscCall(DSSolve(svd->ds,w,NULL));
 61:   PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
 62:   PetscCall(DSSynchronize(svd->ds,w,NULL));

 64:   /* copy singular vectors */
 65:   PetscCall(DSGetArray(svd->ds,DS_MAT_U,&pU));
 66:   PetscCall(DSGetArray(svd->ds,DS_MAT_V,&pV));
 67:   for (i=0;i<n;i++) {
 68:     if (svd->which == SVD_SMALLEST) k = n - i - 1;
 69:     else k = i;
 70:     svd->sigma[k] = PetscRealPart(w[i]);
 71:     PetscCall(BVGetColumn(svd->U,k,&u));
 72:     PetscCall(BVGetColumn(svd->V,k,&v));
 73:     PetscCall(VecGetOwnershipRange(u,&lowu,&highu));
 74:     PetscCall(VecGetOwnershipRange(v,&lowv,&highv));
 75:     PetscCall(VecGetArray(u,&pu));
 76:     PetscCall(VecGetArray(v,&pv));
 77:     if (M>=N) {
 78:       for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
 79:       for (j=lowv;j<highv;j++) pv[j-lowv] = pV[i*ld+j];
 80:     } else {
 81:       for (j=lowu;j<highu;j++) pu[j-lowu] = pV[i*ld+j];
 82:       for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
 83:     }
 84:     PetscCall(VecRestoreArray(u,&pu));
 85:     PetscCall(VecRestoreArray(v,&pv));
 86:     PetscCall(BVRestoreColumn(svd->U,k,&u));
 87:     PetscCall(BVRestoreColumn(svd->V,k,&v));
 88:   }
 89:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&pU));
 90:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&pV));

 92:   svd->nconv  = n;
 93:   svd->its    = 1;
 94:   svd->reason = SVD_CONVERGED_TOL;

 96:   PetscCall(MatDestroy(&mat));
 97:   PetscCall(PetscFree(w));
 98:   PetscFunctionReturn(PETSC_SUCCESS);
 99: }

101: static PetscErrorCode SVDSolve_LAPACK_GSVD(SVD svd)
102: {
103:   PetscInt          nsv,m,n,p,i,j,mlocal,plocal,ld,lowx,lowu,lowv,highx;
104:   Mat               Ar,A,Ads,Br,B,Bds;
105:   Vec               uv,x;
106:   PetscScalar       *U,*V,*X,*px,*puv,*w;

108:   PetscFunctionBegin;
109:   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
110:   PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
111:   PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&A));
112:   PetscCall(MatDestroy(&Ar));
113:   PetscCall(MatCreateRedundantMatrix(svd->OPb,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Br));
114:   PetscCall(MatConvert(Br,MATSEQDENSE,MAT_INITIAL_MATRIX,&B));
115:   PetscCall(MatDestroy(&Br));
116:   PetscCall(MatGetSize(A,&m,&n));
117:   PetscCall(MatGetLocalSize(svd->OP,&mlocal,NULL));
118:   PetscCall(MatGetLocalSize(svd->OPb,&plocal,NULL));
119:   PetscCall(MatGetSize(B,&p,NULL));
120:   PetscCall(DSSetDimensions(svd->ds,m,0,0));
121:   PetscCall(DSGSVDSetDimensions(svd->ds,n,p));
122:   PetscCall(DSGetMat(svd->ds,DS_MAT_A,&Ads));
123:   PetscCall(MatCopy(A,Ads,SAME_NONZERO_PATTERN));
124:   PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&Ads));
125:   PetscCall(DSGetMat(svd->ds,DS_MAT_B,&Bds));
126:   PetscCall(MatCopy(B,Bds,SAME_NONZERO_PATTERN));
127:   PetscCall(DSRestoreMat(svd->ds,DS_MAT_B,&Bds));
128:   PetscCall(DSSetState(svd->ds,DS_STATE_RAW));

130:   nsv = PetscMin(n,PetscMin(p,m));
131:   PetscCall(PetscMalloc1(nsv,&w));
132:   PetscCall(DSSolve(svd->ds,w,NULL));
133:   PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
134:   PetscCall(DSSynchronize(svd->ds,w,NULL));
135:   PetscCall(DSGetDimensions(svd->ds,NULL,NULL,NULL,&nsv));

137:   /* copy singular vectors */
138:   PetscCall(MatGetOwnershipRange(svd->OP,&lowu,NULL));
139:   PetscCall(MatGetOwnershipRange(svd->OPb,&lowv,NULL));
140:   PetscCall(DSGetArray(svd->ds,DS_MAT_X,&X));
141:   PetscCall(DSGetArray(svd->ds,DS_MAT_U,&U));
142:   PetscCall(DSGetArray(svd->ds,DS_MAT_V,&V));
143:   for (j=0;j<nsv;j++) {
144:     svd->sigma[j] = PetscRealPart(w[j]);
145:     PetscCall(BVGetColumn(svd->V,j,&x));
146:     PetscCall(VecGetOwnershipRange(x,&lowx,&highx));
147:     PetscCall(VecGetArrayWrite(x,&px));
148:     for (i=lowx;i<highx;i++) px[i-lowx] = X[i+j*ld];
149:     PetscCall(VecRestoreArrayWrite(x,&px));
150:     PetscCall(BVRestoreColumn(svd->V,j,&x));
151:     PetscCall(BVGetColumn(svd->U,j,&uv));
152:     PetscCall(VecGetArrayWrite(uv,&puv));
153:     for (i=0;i<mlocal;i++) puv[i] = U[i+lowu+j*ld];
154:     for (i=0;i<plocal;i++) puv[i+mlocal] = V[i+lowv+j*ld];
155:     PetscCall(VecRestoreArrayWrite(uv,&puv));
156:     PetscCall(BVRestoreColumn(svd->U,j,&uv));
157:   }
158:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_X,&X));
159:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&U));
160:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&V));

162:   svd->nconv  = nsv;
163:   svd->its    = 1;
164:   svd->reason = SVD_CONVERGED_TOL;

166:   PetscCall(MatDestroy(&A));
167:   PetscCall(MatDestroy(&B));
168:   PetscCall(PetscFree(w));
169:   PetscFunctionReturn(PETSC_SUCCESS);
170: }

172: static PetscErrorCode SVDSolve_LAPACK_HSVD(SVD svd)
173: {
174:   PetscInt          M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
175:   Mat               A,Ar,mat,D;
176:   Vec               u,v,vomega;
177:   PetscScalar       *pU,*pV,*pu,*pv,*w;
178:   PetscReal         *pD;

180:   PetscFunctionBegin;
181:   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
182:   PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
183:   PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat));
184:   PetscCall(MatDestroy(&Ar));
185:   PetscCall(MatGetSize(mat,&M,&N));
186:   PetscCall(DSSetDimensions(svd->ds,M,0,0));
187:   PetscCall(DSHSVDSetDimensions(svd->ds,N));
188:   PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
189:   PetscCall(MatCopy(mat,A,SAME_NONZERO_PATTERN));
190:   PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
191:   PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
192:   PetscCall(VecCopy(svd->omega,vomega));
193:   PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
194:   PetscCall(DSSetState(svd->ds,DS_STATE_RAW));

196:   n = PetscMin(M,N);
197:   PetscCall(PetscMalloc1(n,&w));
198:   PetscCall(DSSolve(svd->ds,w,NULL));
199:   PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
200:   PetscCall(DSSynchronize(svd->ds,w,NULL));

202:   /* copy singular vectors */
203:   PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&pD));
204:   PetscCall(DSGetArray(svd->ds,DS_MAT_U,&pU));
205:   PetscCall(DSGetArray(svd->ds,DS_MAT_V,&pV));
206:   for (i=0;i<n;i++) {
207:     if (svd->which == SVD_SMALLEST) k = n - i - 1;
208:     else k = i;
209:     svd->sigma[k] = PetscRealPart(w[i]);
210:     svd->sign[k]  = pD[i];
211:     PetscCall(BVGetColumn(svd->U,k,&u));
212:     PetscCall(BVGetColumn(svd->V,k,&v));
213:     PetscCall(VecGetOwnershipRange(u,&lowu,&highu));
214:     PetscCall(VecGetOwnershipRange(v,&lowv,&highv));
215:     PetscCall(VecGetArray(u,&pu));
216:     PetscCall(VecGetArray(v,&pv));
217:     if (M>=N) {
218:       for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
219:       for (j=lowv;j<highv;j++) pv[j-lowv] = pV[i*ld+j];
220:     } else {
221:       for (j=lowu;j<highu;j++) pu[j-lowu] = pV[i*ld+j];
222:       for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
223:     }
224:     PetscCall(VecRestoreArray(u,&pu));
225:     PetscCall(VecRestoreArray(v,&pv));
226:     PetscCall(BVRestoreColumn(svd->U,k,&u));
227:     PetscCall(BVRestoreColumn(svd->V,k,&v));
228:   }
229:   PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&pD));
230:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&pU));
231:   PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&pV));

233:   svd->nconv  = n;
234:   svd->its    = 1;
235:   svd->reason = SVD_CONVERGED_TOL;

237:   PetscCall(MatDestroy(&mat));
238:   PetscCall(PetscFree(w));
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: static PetscErrorCode SVDSetDSType_LAPACK(SVD svd)
243: {
244:   PetscFunctionBegin;
245:   PetscCall(DSSetType(svd->ds,svd->OPb?DSGSVD:svd->omega?DSHSVD:DSSVD));
246:   PetscFunctionReturn(PETSC_SUCCESS);
247: }

249: SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
250: {
251:   PetscFunctionBegin;
252:   svd->ops->setup     = SVDSetUp_LAPACK;
253:   svd->ops->solve     = SVDSolve_LAPACK;
254:   svd->ops->solveg    = SVDSolve_LAPACK_GSVD;
255:   svd->ops->solveh    = SVDSolve_LAPACK_HSVD;
256:   svd->ops->setdstype = SVDSetDSType_LAPACK;
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }