Actual source code: svdlapack.c
slepc-3.15.2 2021-09-20
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: {
20: PetscInt M,N,P;
23: MatGetSize(svd->A,&M,&N);
24: if (!svd->isgeneralized) svd->ncv = N;
25: else {
26: MatGetSize(svd->OPb,&P,NULL);
27: svd->ncv = PetscMin(M,PetscMin(N,P));
28: }
29: if (svd->mpd!=PETSC_DEFAULT) { PetscInfo(svd,"Warning: parameter mpd ignored\n"); }
30: SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
31: if (svd->max_it==PETSC_DEFAULT) svd->max_it = 1;
32: svd->leftbasis = PETSC_TRUE;
33: SVDAllocateSolution(svd,0);
34: DSSetType(svd->ds,DSSVD);
35: DSAllocate(svd->ds,PetscMax(M,N));
36: return(0);
37: }
39: PetscErrorCode SVDSolve_LAPACK(SVD svd)
40: {
41: PetscErrorCode ierr;
42: PetscInt M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
43: Mat Ar,mat;
44: Vec u,v;
45: PetscScalar *pU,*pVT,*pu,*pv,*A,*w;
46: const PetscScalar *pmat;
49: DSGetLeadingDimension(svd->ds,&ld);
50: MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar);
51: MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat);
52: MatDestroy(&Ar);
53: MatGetSize(mat,&M,&N);
54: DSSetDimensions(svd->ds,M,N,0,0);
55: MatDenseGetArrayRead(mat,&pmat);
56: DSGetArray(svd->ds,DS_MAT_A,&A);
57: for (i=0;i<M;i++)
58: for (j=0;j<N;j++)
59: A[i+j*ld] = pmat[i+j*M];
60: DSRestoreArray(svd->ds,DS_MAT_A,&A);
61: MatDenseRestoreArrayRead(mat,&pmat);
62: DSSetState(svd->ds,DS_STATE_RAW);
64: n = PetscMin(M,N);
65: PetscMalloc1(n,&w);
66: DSSolve(svd->ds,w,NULL);
67: DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
68: DSSynchronize(svd->ds,w,NULL);
70: /* copy singular vectors */
71: DSGetArray(svd->ds,DS_MAT_U,&pU);
72: DSGetArray(svd->ds,DS_MAT_VT,&pVT);
73: for (i=0;i<n;i++) {
74: if (svd->which == SVD_SMALLEST) k = n - i - 1;
75: else k = i;
76: svd->sigma[k] = PetscRealPart(w[i]);
77: BVGetColumn(svd->U,k,&u);
78: BVGetColumn(svd->V,k,&v);
79: VecGetOwnershipRange(u,&lowu,&highu);
80: VecGetOwnershipRange(v,&lowv,&highv);
81: VecGetArray(u,&pu);
82: VecGetArray(v,&pv);
83: if (M>=N) {
84: for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
85: for (j=lowv;j<highv;j++) pv[j-lowv] = PetscConj(pVT[j*ld+i]);
86: } else {
87: for (j=lowu;j<highu;j++) pu[j-lowu] = PetscConj(pVT[j*ld+i]);
88: for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
89: }
90: VecRestoreArray(u,&pu);
91: VecRestoreArray(v,&pv);
92: BVRestoreColumn(svd->U,k,&u);
93: BVRestoreColumn(svd->V,k,&v);
94: }
95: DSRestoreArray(svd->ds,DS_MAT_U,&pU);
96: DSRestoreArray(svd->ds,DS_MAT_VT,&pVT);
98: svd->nconv = n;
99: svd->reason = SVD_CONVERGED_TOL;
101: MatDestroy(&mat);
102: PetscFree(w);
103: return(0);
104: }
106: PetscErrorCode SVDSolve_LAPACK_GSVD(SVD svd)
107: {
109: PetscInt m,n,p,i,j,mlocal,plocal,ld,lowx,lowu,lowv,highx;
110: PetscBLASInt m_,n_,p_,q_,r_,k,l,lda_,ldb_,ldu_,ldv_,ldq_,lwork,*iwork,info;
111: Mat Ar,A,Br,B;
112: Vec uv,x;
113: PetscScalar *pA,*pB,*U,*V,*Q,*px,*puv,*work,sone=1.0,smone=-1.0;
114: PetscReal *alpha,*beta;
115: #if defined (PETSC_USE_COMPLEX)
116: PetscReal *rwork;
117: #endif
118: #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
119: PetscScalar a,dummy;
120: PetscReal rdummy;
121: PetscBLASInt idummy;
122: #endif
125: DSGetLeadingDimension(svd->ds,&ld);
126: MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar);
127: MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&A);
128: MatDestroy(&Ar);
129: MatCreateRedundantMatrix(svd->OPb,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Br);
130: MatConvert(Br,MATSEQDENSE,MAT_INITIAL_MATRIX,&B);
131: MatDestroy(&Br);
132: MatGetSize(A,&m,&n);
133: MatGetLocalSize(svd->OP,&mlocal,NULL);
134: MatGetLocalSize(svd->OPb,&plocal,NULL);
135: MatGetSize(B,&p,NULL);
136: PetscBLASIntCast(m,&m_);
137: PetscBLASIntCast(n,&n_);
138: PetscBLASIntCast(p,&p_);
139: lda_ = m_; ldb_ = p_;
140: ldu_ = m_; ldv_ = p_; ldq_ = n_;
142: MatDenseGetArray(A,&pA);
143: MatDenseGetArray(B,&pB);
145: #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
147: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
148: /* workspace query and memory allocation */
149: lwork = -1;
150: #if !defined (PETSC_USE_COMPLEX)
151: PetscStackCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m_,&n_,&p_,&k,&l,&dummy,&lda_,&dummy,&ldb_,&rdummy,&rdummy,&dummy,&ldu_,&dummy,&ldv_,&dummy,&ldq_,&a,&lwork,&idummy,&info));
152: PetscBLASIntCast((PetscInt)a,&lwork);
153: #else
154: PetscStackCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m_,&n_,&p_,&k,&l,&dummy,&lda_,&dummy,&ldb_,&rdummy,&rdummy,&dummy,&ldu_,&dummy,&ldv_,&dummy,&ldq_,&a,&lwork,&rdummy,&idummy,&info));
155: PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
156: #endif
157: PetscMalloc7(m*m,&U,p*p,&V,n*n,&Q,n,&alpha,n,&beta,lwork,&work,n,&iwork);
159: #if !defined (PETSC_USE_COMPLEX)
160: PetscStackCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m_,&n_,&p_,&k,&l,pA,&lda_,pB,&ldb_,alpha,beta,U,&ldu_,V,&ldv_,Q,&ldq_,work,&lwork,iwork,&info));
161: #else
162: PetscMalloc1(2*n,&rwork);
163: PetscStackCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m_,&n_,&p_,&k,&l,pA,&lda_,pB,&ldb_,alpha,beta,U,&ldu_,V,&ldv_,Q,&ldq_,work,&lwork,rwork,iwork,&info));
164: #endif
165: PetscFPTrapPop();
166: SlepcCheckLapackInfo("ggsvd3",info);
168: #else // defined(SLEPC_MISSING_LAPACK_GGSVD3)
170: lwork = PetscMax(PetscMax(3*n,m),p)+n;
171: PetscMalloc7(m*m,&U,p*p,&V,n*n,&Q,n,&alpha,n,&beta,lwork,&work,n,&iwork);
173: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
174: #if !defined (PETSC_USE_COMPLEX)
175: PetscStackCallBLAS("LAPACKggsvd",LAPACKggsvd_("U","V","Q",&m_,&n_,&p_,&k,&l,pA,&lda_,pB,&ldb_,alpha,beta,U,&ldu_,V,&ldv_,Q,&ldq_,work,iwork,&info));
176: #else
177: PetscMalloc1(2*n,&rwork);
178: PetscStackCallBLAS("LAPACKggsvd",LAPACKggsvd_("U","V","Q",&m_,&n_,&p_,&k,&l,pA,&lda_,pB,&ldb_,alpha,beta,U,&ldu_,V,&ldv_,Q,&ldq_,work,rwork,iwork,&info));
179: #endif
180: PetscFPTrapPop();
181: SlepcCheckLapackInfo("ggsvd",info);
183: #endif
185: if (k+l<n) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"The case k+l<n not supported yet");
187: /* X = Q*inv(R) */
188: q_ = PetscMin(m_,n_);
189: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n_,&q_,&sone,pA,&lda_,Q,&ldq_));
190: if (m<n) {
191: r_ = n_-m_;
192: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&r_,&m_,&sone,Q,&ldq_,pA,&lda_,&smone,Q+m_*ldq_,&ldq_));
193: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n_,&r_,&sone,pB+m_*ldb_,&ldb_,Q+m_*ldq_,&ldq_));
194: }
195: MatDenseRestoreArray(A,&pA);
196: MatDenseRestoreArray(B,&pB);
198: /* copy singular triplets */
199: for (j=k;j<PetscMin(m,k+l);j++) {
200: svd->sigma[j-k] = alpha[j]/beta[j];
201: BVGetColumn(svd->V,j-k,&x);
202: VecGetOwnershipRange(x,&lowx,&highx);
203: VecGetArrayWrite(x,&px);
204: for (i=lowx;i<highx;i++) px[i-lowx] = Q[i+j*n];
205: VecRestoreArrayWrite(x,&px);
206: BVRestoreColumn(svd->V,j-k,&x);
207: BVGetColumn(svd->U,j-k,&uv);
208: MatGetOwnershipRange(svd->OP,&lowu,NULL);
209: MatGetOwnershipRange(svd->OPb,&lowv,NULL);
210: VecGetArrayWrite(uv,&puv);
211: for (i=0;i<mlocal;i++) puv[i] = U[i+lowu+j*m];
212: for (i=0;i<plocal;i++) puv[i+mlocal] = V[i+lowv+(j-k)*p];
213: VecRestoreArrayWrite(uv,&puv);
214: BVRestoreColumn(svd->U,j-k,&uv);
215: }
217: svd->nconv = PetscMin(m,k+l)-k;
218: svd->reason = SVD_CONVERGED_TOL;
220: MatDestroy(&A);
221: MatDestroy(&B);
222: PetscFree7(U,V,Q,alpha,beta,work,iwork);
223: #if defined (PETSC_USE_COMPLEX)
224: PetscFree(rwork);
225: #endif
226: return(0);
227: }
229: SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
230: {
232: svd->ops->setup = SVDSetUp_LAPACK;
233: svd->ops->solve = SVDSolve_LAPACK;
234: svd->ops->solveg = SVDSolve_LAPACK_GSVD;
235: return(0);
236: }