Actual source code: svdlapack.c

slepc-3.18.2 2023-01-26
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: PetscErrorCode SVDSetUp_LAPACK(SVD svd)
 18: {
 19:   PetscInt       M,N,P=0;

 21:   SVDCheckDefinite(svd);
 22:   MatGetSize(svd->A,&M,&N);
 23:   if (!svd->isgeneralized) svd->ncv = N;
 24:   else {
 25:     MatGetSize(svd->OPb,&P,NULL);
 26:     svd->ncv = PetscMin(M,PetscMin(N,P));
 27:   }
 28:   if (svd->mpd!=PETSC_DEFAULT) PetscInfo(svd,"Warning: parameter mpd ignored\n");
 29:   SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
 30:   if (svd->max_it==PETSC_DEFAULT) svd->max_it = 1;
 31:   svd->leftbasis = PETSC_TRUE;
 32:   SVDAllocateSolution(svd,0);
 33:   DSSetType(svd->ds,svd->isgeneralized?DSGSVD:DSSVD);
 34:   DSAllocate(svd->ds,PetscMax(N,PetscMax(M,P)));
 35:   return 0;
 36: }

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

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

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

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

 95:   MatDestroy(&mat);
 96:   PetscFree(w);
 97:   return 0;
 98: }

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

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

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

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

160:   svd->nconv  = nsv;
161:   svd->its    = 1;
162:   svd->reason = SVD_CONVERGED_TOL;

164:   MatDestroy(&A);
165:   MatDestroy(&B);
166:   PetscFree(w);
167:   return 0;
168: }

170: SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
171: {
172:   svd->ops->setup   = SVDSetUp_LAPACK;
173:   svd->ops->solve   = SVDSolve_LAPACK;
174:   svd->ops->solveg  = SVDSolve_LAPACK_GSVD;
175:   return 0;
176: }