@@ -50,16 +50,26 @@ object EigenValueDecomposition {
50
50
51
51
val arpack = ARPACK .getInstance()
52
52
53
+ // tolerance used in stopping criterion
53
54
val tolW = new doubleW(tol)
55
+ // number of desired eigenvalues, 0 < nev < n
54
56
val nev = new intW(k)
57
+ // nev Lanczos vectors are generated are generated in the first iteration
58
+ // ncv-nev Lanczos vectors are generated are generated in each subsequent iteration
59
+ // ncv must be smaller than n
55
60
val ncv = scala.math.min(2 * k, n)
56
61
62
+ // "I" for standard eigenvalue problem, "G" for generalized eigenvalue problem
57
63
val bmat = " I"
64
+ // "LM" : compute the NEV largest (in magnitude) eigenvalues
58
65
val which = " LM"
59
66
60
67
var iparam = new Array [Int ](11 )
68
+ // use exact shift in each iteration
61
69
iparam(0 ) = 1
70
+ // maximum number of Arnoldi update iterations, or the actual number of iterations on output
62
71
iparam(2 ) = 300
72
+ // Mode 1: A*x = lambda*x, A symmetric
63
73
iparam(6 ) = 1
64
74
65
75
var ido = new intW(0 )
@@ -70,44 +80,58 @@ object EigenValueDecomposition {
70
80
var workl = new Array [Double ](ncv * (ncv + 8 ))
71
81
var ipntr = new Array [Int ](11 )
72
82
73
- // first call to ARPACK
83
+ // call ARPACK's reverse communication, first iteration with ido = 0
74
84
arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr, workd,
75
85
workl, workl.length, info)
76
86
77
87
val w = BDV (workd)
78
88
79
- while (ido.`val` != 99 ) {
89
+ // ido = 99 : done flag in reverse communication
90
+ while (ido.`val` != 99 ) {
80
91
if (ido.`val` != - 1 && ido.`val` != 1 ) {
81
- throw new IllegalStateException (" ARPACK returns ido = " + ido.`val`)
92
+ throw new IllegalStateException (" ARPACK returns ido = " + ido.`val` +
93
+ " This flag is not compatible with Mode 1: A*x = lambda*x, A symmetric." )
82
94
}
83
95
// multiply working vector with the matrix
84
96
val inputOffset = ipntr(0 ) - 1
85
97
val outputOffset = ipntr(1 ) - 1
86
98
val x = w(inputOffset until inputOffset + n)
87
99
val y = w(outputOffset until outputOffset + n)
88
100
y := BDV (mul(Vectors .fromBreeze(x).asInstanceOf [DenseVector ]).toArray)
89
- // call ARPACK
101
+ // call ARPACK's reverse communication
90
102
arpack.dsaupd(ido, bmat, n, which, nev.`val`, tolW, resid, ncv, v, n, iparam, ipntr,
91
103
workd, workl, workl.length, info)
92
104
}
93
105
94
106
if (info.`val` != 0 ) {
95
- throw new IllegalStateException (" ARPACK returns non-zero info = " + info.`val`)
107
+ info.`val` match {
108
+ case 1 => throw new IllegalStateException (" ARPACK returns non-zero info = " + info.`val` +
109
+ " Maximum number of iterations taken. (Refer ARPACK user guide for details)" )
110
+ case 2 => throw new IllegalStateException (" ARPACK returns non-zero info = " + info.`val` +
111
+ " No shifts could be applied. Try to increase NCV. " +
112
+ " (Refer ARPACK user guide for details)" )
113
+ case _ => throw new IllegalStateException (" ARPACK returns non-zero info = " + info.`val` +
114
+ " Please refer ARPACK user guide for error message." )
115
+ }
96
116
}
97
117
98
118
val d = new Array [Double ](nev.`val`)
99
119
val select = new Array [Boolean ](ncv)
120
+ // copy the Ritz vectors
100
121
val z = java.util.Arrays .copyOfRange(v, 0 , nev.`val` * n)
101
122
123
+ // call ARPACK's post-processing for eigenvectors
102
124
arpack.dseupd(true , " A" , select, d, z, n, 0.0 , bmat, n, which, nev, tol, resid, ncv, v, n,
103
125
iparam, ipntr, workd, workl, workl.length, info)
104
126
127
+ // number of computed eigenvalues, might be smaller than k
105
128
val computed = iparam(4 )
106
129
107
130
val eigenPairs = java.util.Arrays .copyOfRange(d, 0 , computed).zipWithIndex.map{
108
131
r => (r._1, java.util.Arrays .copyOfRange(z, r._2 * n, r._2 * n + n))
109
132
}
110
133
134
+ // sort the eigen-pairs in descending order
111
135
val sortedEigenPairs = eigenPairs.sortBy(- 1 * _._1)
112
136
113
137
// copy eigenvectors in descending order of eigenvalues
0 commit comments