Actual source code: svdlapack.c
slepc-3.21.1 2024-04-26
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: PetscCall(MatGetSize(svd->A,&M,&N));
23: if (!svd->isgeneralized) svd->ncv = N;
24: else {
25: PetscCall(MatGetSize(svd->OPb,&P,NULL));
26: svd->ncv = PetscMin(M,PetscMin(N,P));
27: }
28: if (svd->mpd!=PETSC_DEFAULT) PetscCall(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: PetscCall(SVDAllocateSolution(svd,0));
33: PetscCall(DSAllocate(svd->ds,PetscMax(N,PetscMax(M,P))));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode SVDSolve_LAPACK(SVD svd)
38: {
39: PetscInt M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
40: Mat A,Ar,mat;
41: Vec u,v;
42: PetscScalar *pU,*pV,*pu,*pv,*w;
44: PetscFunctionBegin;
45: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
46: PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
47: PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat));
48: PetscCall(MatDestroy(&Ar));
49: PetscCall(MatGetSize(mat,&M,&N));
50: PetscCall(DSSetDimensions(svd->ds,M,0,0));
51: PetscCall(DSSVDSetDimensions(svd->ds,N));
52: PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
53: PetscCall(MatCopy(mat,A,SAME_NONZERO_PATTERN));
54: PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
55: PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
57: n = PetscMin(M,N);
58: PetscCall(PetscMalloc1(n,&w));
59: PetscCall(DSSolve(svd->ds,w,NULL));
60: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
61: PetscCall(DSSynchronize(svd->ds,w,NULL));
63: /* copy singular vectors */
64: PetscCall(DSGetArray(svd->ds,DS_MAT_U,&pU));
65: PetscCall(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: PetscCall(BVGetColumn(svd->U,k,&u));
71: PetscCall(BVGetColumn(svd->V,k,&v));
72: PetscCall(VecGetOwnershipRange(u,&lowu,&highu));
73: PetscCall(VecGetOwnershipRange(v,&lowv,&highv));
74: PetscCall(VecGetArray(u,&pu));
75: PetscCall(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: PetscCall(VecRestoreArray(u,&pu));
84: PetscCall(VecRestoreArray(v,&pv));
85: PetscCall(BVRestoreColumn(svd->U,k,&u));
86: PetscCall(BVRestoreColumn(svd->V,k,&v));
87: }
88: PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&pU));
89: PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&pV));
91: svd->nconv = n;
92: svd->its = 1;
93: svd->reason = SVD_CONVERGED_TOL;
95: PetscCall(MatDestroy(&mat));
96: PetscCall(PetscFree(w));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static 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: PetscFunctionBegin;
108: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
109: PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
110: PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&A));
111: PetscCall(MatDestroy(&Ar));
112: PetscCall(MatCreateRedundantMatrix(svd->OPb,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Br));
113: PetscCall(MatConvert(Br,MATSEQDENSE,MAT_INITIAL_MATRIX,&B));
114: PetscCall(MatDestroy(&Br));
115: PetscCall(MatGetSize(A,&m,&n));
116: PetscCall(MatGetLocalSize(svd->OP,&mlocal,NULL));
117: PetscCall(MatGetLocalSize(svd->OPb,&plocal,NULL));
118: PetscCall(MatGetSize(B,&p,NULL));
119: PetscCall(DSSetDimensions(svd->ds,m,0,0));
120: PetscCall(DSGSVDSetDimensions(svd->ds,n,p));
121: PetscCall(DSGetMat(svd->ds,DS_MAT_A,&Ads));
122: PetscCall(MatCopy(A,Ads,SAME_NONZERO_PATTERN));
123: PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&Ads));
124: PetscCall(DSGetMat(svd->ds,DS_MAT_B,&Bds));
125: PetscCall(MatCopy(B,Bds,SAME_NONZERO_PATTERN));
126: PetscCall(DSRestoreMat(svd->ds,DS_MAT_B,&Bds));
127: PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
129: nsv = PetscMin(n,PetscMin(p,m));
130: PetscCall(PetscMalloc1(nsv,&w));
131: PetscCall(DSSolve(svd->ds,w,NULL));
132: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
133: PetscCall(DSSynchronize(svd->ds,w,NULL));
134: PetscCall(DSGetDimensions(svd->ds,NULL,NULL,NULL,&nsv));
136: /* copy singular vectors */
137: PetscCall(MatGetOwnershipRange(svd->OP,&lowu,NULL));
138: PetscCall(MatGetOwnershipRange(svd->OPb,&lowv,NULL));
139: PetscCall(DSGetArray(svd->ds,DS_MAT_X,&X));
140: PetscCall(DSGetArray(svd->ds,DS_MAT_U,&U));
141: PetscCall(DSGetArray(svd->ds,DS_MAT_V,&V));
142: for (j=0;j<nsv;j++) {
143: svd->sigma[j] = PetscRealPart(w[j]);
144: PetscCall(BVGetColumn(svd->V,j,&x));
145: PetscCall(VecGetOwnershipRange(x,&lowx,&highx));
146: PetscCall(VecGetArrayWrite(x,&px));
147: for (i=lowx;i<highx;i++) px[i-lowx] = X[i+j*ld];
148: PetscCall(VecRestoreArrayWrite(x,&px));
149: PetscCall(BVRestoreColumn(svd->V,j,&x));
150: PetscCall(BVGetColumn(svd->U,j,&uv));
151: PetscCall(VecGetArrayWrite(uv,&puv));
152: for (i=0;i<mlocal;i++) puv[i] = U[i+lowu+j*ld];
153: for (i=0;i<plocal;i++) puv[i+mlocal] = V[i+lowv+j*ld];
154: PetscCall(VecRestoreArrayWrite(uv,&puv));
155: PetscCall(BVRestoreColumn(svd->U,j,&uv));
156: }
157: PetscCall(DSRestoreArray(svd->ds,DS_MAT_X,&X));
158: PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&U));
159: PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&V));
161: svd->nconv = nsv;
162: svd->its = 1;
163: svd->reason = SVD_CONVERGED_TOL;
165: PetscCall(MatDestroy(&A));
166: PetscCall(MatDestroy(&B));
167: PetscCall(PetscFree(w));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode SVDSolve_LAPACK_HSVD(SVD svd)
172: {
173: PetscInt M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
174: Mat A,Ar,mat,D;
175: Vec u,v,vomega;
176: PetscScalar *pU,*pV,*pu,*pv,*w;
177: PetscReal *pD;
179: PetscFunctionBegin;
180: PetscCall(DSGetLeadingDimension(svd->ds,&ld));
181: PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
182: PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat));
183: PetscCall(MatDestroy(&Ar));
184: PetscCall(MatGetSize(mat,&M,&N));
185: PetscCall(DSSetDimensions(svd->ds,M,0,0));
186: PetscCall(DSHSVDSetDimensions(svd->ds,N));
187: PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
188: PetscCall(MatCopy(mat,A,SAME_NONZERO_PATTERN));
189: PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
190: PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
191: PetscCall(VecCopy(svd->omega,vomega));
192: PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
193: PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
195: n = PetscMin(M,N);
196: PetscCall(PetscMalloc1(n,&w));
197: PetscCall(DSSolve(svd->ds,w,NULL));
198: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
199: PetscCall(DSSynchronize(svd->ds,w,NULL));
201: /* copy singular vectors */
202: PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&pD));
203: PetscCall(DSGetArray(svd->ds,DS_MAT_U,&pU));
204: PetscCall(DSGetArray(svd->ds,DS_MAT_V,&pV));
205: for (i=0;i<n;i++) {
206: if (svd->which == SVD_SMALLEST) k = n - i - 1;
207: else k = i;
208: svd->sigma[k] = PetscRealPart(w[i]);
209: svd->sign[k] = pD[i];
210: PetscCall(BVGetColumn(svd->U,k,&u));
211: PetscCall(BVGetColumn(svd->V,k,&v));
212: PetscCall(VecGetOwnershipRange(u,&lowu,&highu));
213: PetscCall(VecGetOwnershipRange(v,&lowv,&highv));
214: PetscCall(VecGetArray(u,&pu));
215: PetscCall(VecGetArray(v,&pv));
216: if (M>=N) {
217: for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
218: for (j=lowv;j<highv;j++) pv[j-lowv] = pV[i*ld+j];
219: } else {
220: for (j=lowu;j<highu;j++) pu[j-lowu] = pV[i*ld+j];
221: for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
222: }
223: PetscCall(VecRestoreArray(u,&pu));
224: PetscCall(VecRestoreArray(v,&pv));
225: PetscCall(BVRestoreColumn(svd->U,k,&u));
226: PetscCall(BVRestoreColumn(svd->V,k,&v));
227: }
228: PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&pD));
229: PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&pU));
230: PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&pV));
232: svd->nconv = n;
233: svd->its = 1;
234: svd->reason = SVD_CONVERGED_TOL;
236: PetscCall(MatDestroy(&mat));
237: PetscCall(PetscFree(w));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode SVDSetDSType_LAPACK(SVD svd)
242: {
243: PetscFunctionBegin;
244: PetscCall(DSSetType(svd->ds,svd->OPb?DSGSVD:svd->omega?DSHSVD:DSSVD));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
249: {
250: PetscFunctionBegin;
251: svd->ops->setup = SVDSetUp_LAPACK;
252: svd->ops->solve = SVDSolve_LAPACK;
253: svd->ops->solveg = SVDSolve_LAPACK_GSVD;
254: svd->ops->solveh = SVDSolve_LAPACK_HSVD;
255: svd->ops->setdstype = SVDSetDSType_LAPACK;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }