1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
/*
ARPACK++ v1.2 2/20/2000
c++ interface to ARPACK code.
MODULE LSVDSol.h
Template functions that exemplify how to print information
about the singular value decomposition of a real nonsymmetric
matrix.
ARPACK Authors
Richard Lehoucq
Danny Sorensen
Chao Yang
Dept. of Computational & Applied Mathematics
Rice University
Houston, Texas
*/
#ifndef LSVDSOL_H
#define LSVDSOL_H
#include <math.h>
#include "blas1c.h"
#include "lapackc.h"
#include "arssym.h"
template<class ARMATRIX, class ARFLOAT, class ARFOP>
void Solution(ARMATRIX &A, ARSymStdEig<ARFLOAT, ARFOP> &Prob)
/*
Prints singular values and vectors on standard "cout" stream.
*/
{
int i, m, n, nAx, nconv;
ARFLOAT *Ax;
ARFLOAT *ResNorm;
m = A.nrows();
n = A.ncols();
nAx = (n>m)?n:m;
nconv = Prob.ConvergedEigenvalues();
cout << endl << endl << "Testing ARPACK++ class ARSymStdEig \n";
cout << "SVD problems: A = U*S*V'" << endl;
cout << "Dimension of the system : " << Prob.GetN() << endl;
cout << "Number of 'requested' eigenvalues : " << Prob.GetNev() << endl;
cout << "Number of 'converged' eigenvalues : " << nconv << endl;
cout << "Number of Arnoldi vectors generated: " << Prob.GetNcv() << endl;
cout << endl;
if (Prob.EigenvaluesFound()) {
// Printing singular values.
cout << "Singular values:" << endl;
for (i=0; i<nconv; i++) {
cout << " sigma[" << (i+1) << "]: " << Prob.Eigenvalue(i) << endl;
}
cout << endl;
}
if (Prob.EigenvectorsFound()) {
// Printing the residual norm || A*v - sigma*u ||
// for the nconv accurately computed vectors u and v.
Ax = new ARFLOAT[nAx];
ResNorm = new ARFLOAT[nconv+1];
// Printing the residual norm || A*v - sigma*u ||
// for the nconv accurately computed vectors u and v.
for (i=0; i<nconv; i++) {
A.MultMv(Prob.RawEigenvector(i)+m, Ax);
axpy(m, -Prob.Eigenvalue(i), Prob.RawEigenvector(i), 1, Ax, 1);
ResNorm[i] = nrm2(m, Ax, 1)/fabs(Prob.Eigenvalue(i));
}
for (i=0; i<nconv; i++) {
cout << "||A*v(" << (i+1) << ") - sigma(" << (i+1);
cout << ")*u(" << (i+1) << ")||: " << ResNorm[i] << endl;
}
cout << endl;
// Printing the residual norm || A'*u - sigma*v ||
// for the nconv accurately computed vectors u and v.
for (i=0; i<nconv; i++) {
A.MultMtv(Prob.RawEigenvector(i), Ax);
axpy(n, -Prob.Eigenvalue(i), Prob.RawEigenvector(i)+m, 1, Ax, 1);
ResNorm[i] = nrm2(n, Ax, 1)/fabs(Prob.Eigenvalue(i));
}
for (i=0; i<nconv; i++) {
cout << "||A'*u(" << (i+1) << ") - sigma(" << (i+1);
cout << ")*v(" << (i+1) << ")||: " << ResNorm[i] << endl;
}
cout << endl;
delete[] Ax;
delete[] ResNorm;
}
} // Solution
#endif // LSVDSOL_H
|
