| 1298 |
slepc |
1 |
/*
|
|
|
2 |
This file contains routines for handling small-size dense problems.
|
|
|
3 |
All routines are simply wrappers to LAPACK routines.
|
| 1376 |
slepc |
4 |
|
|
|
5 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
6 |
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
|
|
7 |
Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain
|
|
|
8 |
|
|
|
9 |
This file is part of SLEPc. See the README file for conditions of use
|
|
|
10 |
and additional information.
|
|
|
11 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
| 1298 |
slepc |
12 |
*/
|
|
|
13 |
|
| 1339 |
slepc |
14 |
#include "src/svd/svdimpl.h" /*I "slepcsvd.h" I*/
|
| 1298 |
slepc |
15 |
#include "slepcblaslapack.h"
|
|
|
16 |
|
| 1333 |
slepc |
17 |
#undef __FUNCT__
|
|
|
18 |
#define __FUNCT__ "SVDDense"
|
|
|
19 |
/*@
|
|
|
20 |
SVDDense - Solves a dense singular value problem.
|
|
|
21 |
|
|
|
22 |
Not Collective
|
|
|
23 |
|
|
|
24 |
Input Parameters:
|
|
|
25 |
+ M - dimension of the problem (rows)
|
|
|
26 |
. N - dimension of the problem (colums)
|
|
|
27 |
- A - pointer to the array containing the matrix values
|
|
|
28 |
|
|
|
29 |
Output Parameters:
|
|
|
30 |
+ sigma - pointer to the array to store the computed singular values
|
|
|
31 |
. U - pointer to the array to store left singular vectors
|
|
|
32 |
- VT - pointer to the array to store right singular vectors
|
|
|
33 |
|
|
|
34 |
Matrix A is overwritten.
|
|
|
35 |
|
|
|
36 |
This routine uses LAPACK routines xGESDD.
|
|
|
37 |
|
|
|
38 |
Level: developer
|
|
|
39 |
|
|
|
40 |
@*/
|
| 1298 |
slepc |
41 |
PetscErrorCode SVDDense(int M,int N,PetscScalar* A,PetscReal* sigma,PetscScalar* U,PetscScalar* VT)
|
|
|
42 |
{
|
| 1336 |
slepc |
43 |
#if defined(SLEPC_MISSING_LAPACK_GESDD)
|
|
|
44 |
PetscFunctionBegin;
|
|
|
45 |
SETERRQ(PETSC_ERR_SUP,"GESDD - Lapack routine is unavailable.");
|
|
|
46 |
#else
|
| 1298 |
slepc |
47 |
PetscErrorCode ierr;
|
|
|
48 |
PetscScalar qwork,*work;
|
|
|
49 |
int n,info,lwork,*iwork;
|
|
|
50 |
#if defined(PETSC_USE_COMPLEX)
|
|
|
51 |
PetscReal *rwork;
|
|
|
52 |
#endif
|
|
|
53 |
|
|
|
54 |
PetscFunctionBegin;
|
|
|
55 |
/* workspace query & allocation */
|
| 1339 |
slepc |
56 |
ierr = PetscLogEventBegin(SVD_Dense,0,0,0,0);CHKERRQ(ierr);
|
| 1298 |
slepc |
57 |
n = PetscMin(M,N);
|
|
|
58 |
ierr = PetscMalloc(sizeof(int)*8*n,&iwork);CHKERRQ(ierr);
|
|
|
59 |
lwork = -1;
|
|
|
60 |
#if defined(PETSC_USE_COMPLEX)
|
|
|
61 |
ierr = PetscMalloc(sizeof(PetscReal)*(5*n*n+7*n),&rwork);CHKERRQ(ierr);
|
|
|
62 |
LAPACKgesdd_("O",&M,&N,A,&M,sigma,U,&M,VT,&N,&qwork,&lwork,rwork,iwork,&info,1);
|
|
|
63 |
#else
|
|
|
64 |
LAPACKgesdd_("O",&M,&N,A,&M,sigma,U,&M,VT,&N,&qwork,&lwork,iwork,&info,1);
|
|
|
65 |
#endif
|
|
|
66 |
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGESDD %d",info);
|
| 1371 |
slepc |
67 |
lwork = (int)PetscRealPart(qwork);
|
| 1298 |
slepc |
68 |
ierr = PetscMalloc(sizeof(PetscScalar)*lwork,&work);CHKERRQ(ierr);
|
|
|
69 |
|
|
|
70 |
/* computation */
|
|
|
71 |
#if defined(PETSC_USE_COMPLEX)
|
|
|
72 |
LAPACKgesdd_("O",&M,&N,A,&M,sigma,U,&M,VT,&N,work,&lwork,rwork,iwork,&info,1);
|
|
|
73 |
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
|
|
74 |
#else
|
|
|
75 |
LAPACKgesdd_("O",&M,&N,A,&M,sigma,U,&M,VT,&N,work,&lwork,iwork,&info,1);
|
|
|
76 |
#endif
|
|
|
77 |
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGESDD %d",info);
|
|
|
78 |
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
|
|
79 |
ierr = PetscFree(work);CHKERRQ(ierr);
|
| 1339 |
slepc |
80 |
ierr = PetscLogEventEnd(SVD_Dense,0,0,0,0);CHKERRQ(ierr);
|
| 1298 |
slepc |
81 |
PetscFunctionReturn(0);
|
| 1336 |
slepc |
82 |
#endif
|
| 1326 |
slepc |
83 |
}
|