| 1263 |
slepc |
1 |
/*
|
|
|
2 |
This file implements a wrapper to the LAPACK SVD subroutines.
|
| 1376 |
slepc |
3 |
|
|
|
4 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
5 |
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
|
|
6 |
Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain
|
|
|
7 |
|
|
|
8 |
This file is part of SLEPc. See the README file for conditions of use
|
|
|
9 |
and additional information.
|
|
|
10 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
| 1263 |
slepc |
11 |
*/
|
| 1376 |
slepc |
12 |
|
| 1521 |
slepc |
13 |
#include "private/svdimpl.h"
|
| 1263 |
slepc |
14 |
#include "slepcblaslapack.h"
|
|
|
15 |
|
|
|
16 |
#undef __FUNCT__
|
| 1591 |
slepc |
17 |
#define __FUNCT__ "SVDSetUp_LAPACK"
|
|
|
18 |
PetscErrorCode SVDSetUp_LAPACK(SVD svd)
|
| 1263 |
slepc |
19 |
{
|
|
|
20 |
PetscErrorCode ierr;
|
| 1605 |
slepc |
21 |
PetscInt N,i,nloc;
|
|
|
22 |
PetscScalar *pU;
|
| 1263 |
slepc |
23 |
|
|
|
24 |
PetscFunctionBegin;
|
| 1314 |
slepc |
25 |
ierr = SVDMatGetSize(svd,PETSC_NULL,&N);CHKERRQ(ierr);
|
|
|
26 |
svd->ncv = N;
|
| 1591 |
slepc |
27 |
if (svd->mpd) PetscInfo(svd,"Warning: parameter mpd ignored\n");
|
| 1272 |
slepc |
28 |
svd->max_it = 1;
|
| 1314 |
slepc |
29 |
if (svd->ncv!=svd->n) {
|
|
|
30 |
if (svd->U) {
|
| 1605 |
slepc |
31 |
ierr = VecGetArray(svd->U[0],&pU);CHKERRQ(ierr);
|
| 1314 |
slepc |
32 |
for (i=0;i<svd->n;i++) { ierr = VecDestroy(svd->U[i]); CHKERRQ(ierr); }
|
| 1605 |
slepc |
33 |
ierr = PetscFree(pU);CHKERRQ(ierr);
|
| 1314 |
slepc |
34 |
ierr = PetscFree(svd->U);CHKERRQ(ierr);
|
|
|
35 |
}
|
|
|
36 |
ierr = PetscMalloc(sizeof(Vec)*svd->ncv,&svd->U);CHKERRQ(ierr);
|
| 1605 |
slepc |
37 |
ierr = SVDMatGetLocalSize(svd,&nloc,PETSC_NULL);CHKERRQ(ierr);
|
|
|
38 |
ierr = PetscMalloc(svd->ncv*nloc*sizeof(PetscScalar),&pU);CHKERRQ(ierr);
|
|
|
39 |
for (i=0;i<svd->ncv;i++) {
|
|
|
40 |
ierr = VecCreateMPIWithArray(((PetscObject)svd)->comm,nloc,PETSC_DECIDE,pU+i*nloc,&svd->U[i]);CHKERRQ(ierr);
|
|
|
41 |
}
|
| 1314 |
slepc |
42 |
}
|
| 1263 |
slepc |
43 |
PetscFunctionReturn(0);
|
|
|
44 |
}
|
|
|
45 |
|
|
|
46 |
#undef __FUNCT__
|
|
|
47 |
#define __FUNCT__ "SVDSolve_LAPACK"
|
|
|
48 |
PetscErrorCode SVDSolve_LAPACK(SVD svd)
|
|
|
49 |
{
|
|
|
50 |
PetscErrorCode ierr;
|
| 1504 |
slepc |
51 |
PetscInt M,N,n,i,j,k;
|
| 1263 |
slepc |
52 |
Mat mat;
|
| 1298 |
slepc |
53 |
PetscScalar *pU,*pVT,*pmat,*pu,*pv;
|
| 1289 |
slepc |
54 |
PetscReal *sigma;
|
| 1263 |
slepc |
55 |
|
|
|
56 |
PetscFunctionBegin;
|
| 1318 |
slepc |
57 |
ierr = MatConvert(svd->OP,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat);CHKERRQ(ierr);
|
|
|
58 |
ierr = MatGetArray(mat,&pmat);CHKERRQ(ierr);
|
|
|
59 |
ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
|
|
|
60 |
if (M>=N) {
|
|
|
61 |
n = N;
|
|
|
62 |
pU = PETSC_NULL;
|
|
|
63 |
ierr = PetscMalloc(sizeof(PetscScalar)*N*N,&pVT);CHKERRQ(ierr);
|
| 1263 |
slepc |
64 |
} else {
|
| 1318 |
slepc |
65 |
n = M;
|
|
|
66 |
ierr = PetscMalloc(sizeof(PetscScalar)*M*M,&pU);CHKERRQ(ierr);
|
|
|
67 |
pVT = PETSC_NULL;
|
| 1263 |
slepc |
68 |
}
|
| 1318 |
slepc |
69 |
ierr = PetscMalloc(sizeof(PetscReal)*n,&sigma);CHKERRQ(ierr);
|
| 1263 |
slepc |
70 |
|
| 1298 |
slepc |
71 |
ierr = SVDDense(M,N,pmat,sigma,pU,pVT);CHKERRQ(ierr);
|
|
|
72 |
|
| 1263 |
slepc |
73 |
/* copy singular vectors */
|
| 1318 |
slepc |
74 |
for (i=0;i<n;i++) {
|
|
|
75 |
if (svd->which == SVD_SMALLEST) k = n - i - 1;
|
| 1289 |
slepc |
76 |
else k = i;
|
|
|
77 |
svd->sigma[k] = sigma[i];
|
|
|
78 |
ierr = VecGetArray(svd->U[k],&pu);CHKERRQ(ierr);
|
|
|
79 |
ierr = VecGetArray(svd->V[k],&pv);CHKERRQ(ierr);
|
| 1318 |
slepc |
80 |
if (M>=N) {
|
|
|
81 |
for (j=0;j<M;j++) pu[j] = pmat[i*M+j];
|
|
|
82 |
for (j=0;j<N;j++) pv[j] = pVT[j*N+i];
|
|
|
83 |
} else {
|
|
|
84 |
for (j=0;j<N;j++) pu[j] = pmat[j*M+i];
|
|
|
85 |
for (j=0;j<M;j++) pv[j] = pU[i*M+j];
|
|
|
86 |
}
|
| 1289 |
slepc |
87 |
ierr = VecRestoreArray(svd->U[k],&pu);CHKERRQ(ierr);
|
|
|
88 |
ierr = VecRestoreArray(svd->V[k],&pv);CHKERRQ(ierr);
|
| 1263 |
slepc |
89 |
}
|
|
|
90 |
|
| 1318 |
slepc |
91 |
svd->nconv = n;
|
| 1283 |
slepc |
92 |
svd->reason = SVD_CONVERGED_TOL;
|
|
|
93 |
|
| 1263 |
slepc |
94 |
ierr = MatRestoreArray(mat,&pmat);CHKERRQ(ierr);
|
|
|
95 |
ierr = MatDestroy(mat);CHKERRQ(ierr);
|
| 1289 |
slepc |
96 |
ierr = PetscFree(sigma);CHKERRQ(ierr);
|
| 1318 |
slepc |
97 |
if (M>=N) {
|
|
|
98 |
ierr = PetscFree(pVT);CHKERRQ(ierr);
|
|
|
99 |
} else {
|
|
|
100 |
ierr = PetscFree(pU);CHKERRQ(ierr);
|
|
|
101 |
}
|
| 1263 |
slepc |
102 |
PetscFunctionReturn(0);
|
|
|
103 |
}
|
|
|
104 |
|
|
|
105 |
EXTERN_C_BEGIN
|
|
|
106 |
#undef __FUNCT__
|
|
|
107 |
#define __FUNCT__ "SVDCreate_LAPACK"
|
|
|
108 |
PetscErrorCode SVDCreate_LAPACK(SVD svd)
|
|
|
109 |
{
|
|
|
110 |
PetscFunctionBegin;
|
| 1591 |
slepc |
111 |
svd->ops->setup = SVDSetUp_LAPACK;
|
| 1391 |
slepc |
112 |
svd->ops->solve = SVDSolve_LAPACK;
|
|
|
113 |
svd->ops->destroy = SVDDestroy_Default;
|
| 1283 |
slepc |
114 |
if (svd->transmode == PETSC_DECIDE)
|
|
|
115 |
svd->transmode = SVD_TRANSPOSE_IMPLICIT; /* don't build the transpose */
|
| 1263 |
slepc |
116 |
PetscFunctionReturn(0);
|
|
|
117 |
}
|
|
|
118 |
EXTERN_C_END
|