Actual source code: svdlapack.c
slepc-3.17.2 2022-08-09
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: MatGetSize(svd->A,&M,&N);
22: if (!svd->isgeneralized) svd->ncv = N;
23: else {
24: MatGetSize(svd->OPb,&P,NULL);
25: svd->ncv = PetscMin(M,PetscMin(N,P));
26: }
27: if (svd->mpd!=PETSC_DEFAULT) PetscInfo(svd,"Warning: parameter mpd ignored\n");
28: SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
29: if (svd->max_it==PETSC_DEFAULT) svd->max_it = 1;
30: svd->leftbasis = PETSC_TRUE;
31: SVDAllocateSolution(svd,0);
32: DSSetType(svd->ds,svd->isgeneralized?DSGSVD:DSSVD);
33: DSAllocate(svd->ds,PetscMax(N,PetscMax(M,P)));
34: PetscFunctionReturn(0);
35: }
37: PetscErrorCode SVDSolve_LAPACK(SVD svd)
38: {
39: PetscInt M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
40: Mat Ar,mat;
41: Vec u,v;
42: PetscScalar *pU,*pV,*pu,*pv,*A,*w;
43: const PetscScalar *pmat;
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: MatDenseGetArrayRead(mat,&pmat);
53: DSGetArray(svd->ds,DS_MAT_A,&A);
54: for (i=0;i<M;i++)
55: for (j=0;j<N;j++)
56: A[i+j*ld] = pmat[i+j*M];
57: DSRestoreArray(svd->ds,DS_MAT_A,&A);
58: MatDenseRestoreArrayRead(mat,&pmat);
59: DSSetState(svd->ds,DS_STATE_RAW);
61: n = PetscMin(M,N);
62: PetscMalloc1(n,&w);
63: DSSolve(svd->ds,w,NULL);
64: DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
65: DSSynchronize(svd->ds,w,NULL);
67: /* copy singular vectors */
68: DSGetArray(svd->ds,DS_MAT_U,&pU);
69: DSGetArray(svd->ds,DS_MAT_V,&pV);
70: for (i=0;i<n;i++) {
71: if (svd->which == SVD_SMALLEST) k = n - i - 1;
72: else k = i;
73: svd->sigma[k] = PetscRealPart(w[i]);
74: BVGetColumn(svd->U,k,&u);
75: BVGetColumn(svd->V,k,&v);
76: VecGetOwnershipRange(u,&lowu,&highu);
77: VecGetOwnershipRange(v,&lowv,&highv);
78: VecGetArray(u,&pu);
79: VecGetArray(v,&pv);
80: if (M>=N) {
81: for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
82: for (j=lowv;j<highv;j++) pv[j-lowv] = pV[i*ld+j];
83: } else {
84: for (j=lowu;j<highu;j++) pu[j-lowu] = pV[i*ld+j];
85: for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
86: }
87: VecRestoreArray(u,&pu);
88: VecRestoreArray(v,&pv);
89: BVRestoreColumn(svd->U,k,&u);
90: BVRestoreColumn(svd->V,k,&v);
91: }
92: DSRestoreArray(svd->ds,DS_MAT_U,&pU);
93: DSRestoreArray(svd->ds,DS_MAT_V,&pV);
95: svd->nconv = n;
96: svd->its = 1;
97: svd->reason = SVD_CONVERGED_TOL;
99: MatDestroy(&mat);
100: PetscFree(w);
101: PetscFunctionReturn(0);
102: }
104: PetscErrorCode SVDSolve_LAPACK_GSVD(SVD svd)
105: {
106: PetscInt nsv,m,n,p,i,j,mlocal,plocal,ld,lowx,lowu,lowv,highx;
107: Mat Ar,A,Br,B;
108: Vec uv,x;
109: PetscScalar *Ads,*Bds,*U,*V,*X,*px,*puv,*w;
110: const PetscScalar *pA,*pB;
112: DSGetLeadingDimension(svd->ds,&ld);
113: MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar);
114: MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&A);
115: MatDestroy(&Ar);
116: MatCreateRedundantMatrix(svd->OPb,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Br);
117: MatConvert(Br,MATSEQDENSE,MAT_INITIAL_MATRIX,&B);
118: MatDestroy(&Br);
119: MatGetSize(A,&m,&n);
120: MatGetLocalSize(svd->OP,&mlocal,NULL);
121: MatGetLocalSize(svd->OPb,&plocal,NULL);
122: MatGetSize(B,&p,NULL);
123: DSSetDimensions(svd->ds,m,0,0);
124: DSGSVDSetDimensions(svd->ds,n,p);
125: MatDenseGetArrayRead(A,&pA);
126: MatDenseGetArrayRead(B,&pB);
127: DSGetArray(svd->ds,DS_MAT_A,&Ads);
128: DSGetArray(svd->ds,DS_MAT_B,&Bds);
129: for (j=0;j<n;j++) {
130: for (i=0;i<m;i++) Ads[i+j*ld] = pA[i+j*m];
131: for (i=0;i<p;i++) Bds[i+j*ld] = pB[i+j*p];
132: }
133: DSRestoreArray(svd->ds,DS_MAT_B,&Bds);
134: DSRestoreArray(svd->ds,DS_MAT_A,&Ads);
135: MatDenseRestoreArrayRead(B,&pB);
136: MatDenseRestoreArrayRead(A,&pA);
137: DSSetState(svd->ds,DS_STATE_RAW);
139: nsv = PetscMin(n,PetscMin(p,m));
140: PetscMalloc1(nsv,&w);
141: DSSolve(svd->ds,w,NULL);
142: DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
143: DSSynchronize(svd->ds,w,NULL);
144: DSGetDimensions(svd->ds,NULL,NULL,NULL,&nsv);
146: /* copy singular vectors */
147: MatGetOwnershipRange(svd->OP,&lowu,NULL);
148: MatGetOwnershipRange(svd->OPb,&lowv,NULL);
149: DSGetArray(svd->ds,DS_MAT_X,&X);
150: DSGetArray(svd->ds,DS_MAT_U,&U);
151: DSGetArray(svd->ds,DS_MAT_V,&V);
152: for (j=0;j<nsv;j++) {
153: svd->sigma[j] = PetscRealPart(w[j]);
154: BVGetColumn(svd->V,j,&x);
155: VecGetOwnershipRange(x,&lowx,&highx);
156: VecGetArrayWrite(x,&px);
157: for (i=lowx;i<highx;i++) px[i-lowx] = X[i+j*ld];
158: VecRestoreArrayWrite(x,&px);
159: BVRestoreColumn(svd->V,j,&x);
160: BVGetColumn(svd->U,j,&uv);
161: VecGetArrayWrite(uv,&puv);
162: for (i=0;i<mlocal;i++) puv[i] = U[i+lowu+j*ld];
163: for (i=0;i<plocal;i++) puv[i+mlocal] = V[i+lowv+j*ld];
164: VecRestoreArrayWrite(uv,&puv);
165: BVRestoreColumn(svd->U,j,&uv);
166: }
167: DSRestoreArray(svd->ds,DS_MAT_X,&X);
168: DSRestoreArray(svd->ds,DS_MAT_U,&U);
169: DSRestoreArray(svd->ds,DS_MAT_V,&V);
171: svd->nconv = nsv;
172: svd->its = 1;
173: svd->reason = SVD_CONVERGED_TOL;
175: MatDestroy(&A);
176: MatDestroy(&B);
177: PetscFree(w);
178: PetscFunctionReturn(0);
179: }
181: SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
182: {
183: svd->ops->setup = SVDSetUp_LAPACK;
184: svd->ops->solve = SVDSolve_LAPACK;
185: svd->ops->solveg = SVDSolve_LAPACK_GSVD;
186: PetscFunctionReturn(0);
187: }