Actual source code: svdlapack.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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>
16: PetscErrorCode SVDSetUp_LAPACK(SVD svd)
17: {
19: PetscInt M,N;
22: SVDMatGetSize(svd,&M,&N);
23: svd->ncv = N;
24: if (svd->mpd) { PetscInfo(svd,"Warning: parameter mpd ignored\n"); }
25: if (svd->stop!=SVD_STOP_BASIC) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined stopping test not supported in this solver");
26: svd->max_it = 1;
27: svd->leftbasis = PETSC_TRUE;
28: SVDAllocateSolution(svd,0);
29: DSSetType(svd->ds,DSSVD);
30: DSAllocate(svd->ds,PetscMax(M,N));
31: return(0);
32: }
34: PetscErrorCode SVDSolve_LAPACK(SVD svd)
35: {
37: PetscInt M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
38: Mat Ar,mat;
39: Vec u,v;
40: PetscScalar *pU,*pVT,*pmat,*pu,*pv,*A,*w;
43: DSGetLeadingDimension(svd->ds,&ld);
44: MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar);
45: MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat);
46: MatDestroy(&Ar);
47: MatGetSize(mat,&M,&N);
48: DSSetDimensions(svd->ds,M,N,0,0);
49: MatDenseGetArray(mat,&pmat);
50: DSGetArray(svd->ds,DS_MAT_A,&A);
51: for (i=0;i<M;i++)
52: for (j=0;j<N;j++)
53: A[i+j*ld] = pmat[i+j*M];
54: DSRestoreArray(svd->ds,DS_MAT_A,&A);
55: MatDenseRestoreArray(mat,&pmat);
56: DSSetState(svd->ds,DS_STATE_RAW);
58: n = PetscMin(M,N);
59: PetscMalloc1(n,&w);
60: DSSolve(svd->ds,w,NULL);
61: DSSort(svd->ds,w,NULL,NULL,NULL,NULL);
62: DSSynchronize(svd->ds,w,NULL);
64: /* copy singular vectors */
65: DSGetArray(svd->ds,DS_MAT_U,&pU);
66: DSGetArray(svd->ds,DS_MAT_VT,&pVT);
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: BVGetColumn(svd->U,k,&u);
72: BVGetColumn(svd->V,k,&v);
73: VecGetOwnershipRange(u,&lowu,&highu);
74: VecGetOwnershipRange(v,&lowv,&highv);
75: VecGetArray(u,&pu);
76: 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] = PetscConj(pVT[j*ld+i]);
80: } else {
81: for (j=lowu;j<highu;j++) pu[j-lowu] = PetscConj(pVT[j*ld+i]);
82: for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
83: }
84: VecRestoreArray(u,&pu);
85: VecRestoreArray(v,&pv);
86: BVRestoreColumn(svd->U,k,&u);
87: BVRestoreColumn(svd->V,k,&v);
88: }
89: DSRestoreArray(svd->ds,DS_MAT_U,&pU);
90: DSRestoreArray(svd->ds,DS_MAT_VT,&pVT);
92: svd->nconv = n;
93: svd->reason = SVD_CONVERGED_TOL;
95: MatDestroy(&mat);
96: PetscFree(w);
97: return(0);
98: }
100: SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
101: {
103: svd->ops->setup = SVDSetUp_LAPACK;
104: svd->ops->solve = SVDSolve_LAPACK;
105: return(0);
106: }