Lesson of 24 February 2010¶
Lesson file in a zip file
The header of Jacobi method for eigenvalues computation
File* jacobi.hh
The code of Jacobi method
File jacobi.cc
1#include "jacobi.hh"
2#include "alglin.hh"
3#include <cmath>
4#include <iostream>
5
6namespace Jacobi {
7
8 using namespace std ;
9
10 void
11 GivensLeft( indexType const n, /* Number of rows of matrix A */
12 valueType A[], /* Pointer to coefficients of A */
13 indexType const lda, /* Leading dimension of A, Fortran style */
14 indexType const row1, /* first row for Givens rotation */
15 indexType const row2, /* second row for Givens rotation */
16 valueType const theta ) {
17 //
18 // / \
19 // | cos(theta) sin(theta) |
20 // | |
21 // | -sin(theta) cos(theta) |
22 // \ /
23 //
24 // row1 <- cos(theta) * row1 + sin(theta) * row2
25 // row2 <- -sin(theta) * row1 + cos(theta) * row2
26 valueType cs = cos(theta) ;
27 valueType sn = sin(theta) ;
28 valueType * R1 = A + row1 ;
29 valueType * R2 = A + row2 ;
30 for ( indexType i=0 ; i < n ; ++i, R1 += lda, R2 += lda ) {
31 valueType tmp = cs*(*R1)+sn*(*R2) ;
32 *R2 = cs*(*R2)-sn*(*R1) ;
33 *R1 = tmp ;
34 }
35 }
36
37 void
38 GivensRight( indexType const n, /* Number of rows of matrix A */
39 valueType A[], /* Pointer to coefficients of A */
40 indexType const lda, /* Leading dimension of A, Fortran style */
41 indexType const col1, /* first column for Givens rotation */
42 indexType const col2, /* second column for Givens rotation */
43 valueType const theta ) {
44 //
45 // / \
46 // | cos(theta) -sin(theta) |
47 // | |
48 // | sin(theta) cos(theta) |
49 // \ /
50 //
51 // col1 <- cos(theta) * col1 + sin(theta) * col2
52 // row2 <- -sin(theta) * col1 + cos(theta) * col2
53 valueType cs = cos(theta) ;
54 valueType sn = sin(theta) ;
55 valueType * C1 = A + col1 * lda ;
56 valueType * C2 = A + col2 * lda ;
57 for ( indexType i=0 ; i < n ; ++i, ++C1, ++C2 ) {
58 valueType tmp = cs*(*C1)+sn*(*C2) ;
59 *C2 = cs*(*C2)-sn*(*C1) ;
60 *C1 = tmp ;
61 }
62 }
63
64 valueType
65 evalAngle( indexType const n, /* Number of rows of matrix A */
66 valueType A[], /* Pointer to coefficients of A */
67 indexType const lda, /* Leading dimension of A, Fortran style */
68 indexType const row, // row > col
69 indexType const col ) {
70 //
71 // col
72 // |
73 // . . . . .
74 // . Aii . Aij .
75 // . . . . .
76 // . . . . .
77 // . Aij . Ajj . <--- row
78 // . . . . .
79 valueType Aij = A[row + col * lda] ;
80 valueType Aii = A[col + col * lda] ;
81 //valueType Ajj = A[row + row * lda] ;
82 return -atan2( Aij, Aii ) ;
83 }
84
85 bool
86 Jacobi( indexType const n, /* Number of rows of matrix A */
87 valueType A[], /* Pointer to coefficients of A */
88 indexType const lda, /* Leading dimension of A, Fortran style */
89 valueType Q[], /* Pointer to coefficients of Q */
90 indexType const ldq, /* Leading dimension of Q, Fortran style */
91 valueType const epsi, /* tolerance for off diagonal elements */
92 indexType const maxIter ) { /* maximum number of big iterations */
93 // initialize Q to the identity matrix
94 valueType one = 1 ;
95 alglin::zero( n*ldq, Q, 1 ) ;
96 alglin::copy( n, &one, 0, Q, ldq+1 ) ;
97
98 for ( indexType k = 0 ; k < maxIter ; ++k ) {
99 // evaluate max norm of non diagonal elements
100 indexType iAbs = 0 ;
101 indexType jAbs = 0 ;
102 valueType maxAbs = 0 ;
103 for ( indexType i = 1 ; i < n ; ++i ) {
104 for ( indexType j = 0 ; j < i ; ++j ) {
105 valueType absAij = abs( A[ i + lda * j] ) ;
106 if ( absAij > maxAbs ) {
107 iAbs = i ;
108 jAbs = j ;
109 maxAbs = absAij ;
110 }
111 }
112 }
113 if ( maxAbs < epsi ) return true ;
114 // loop for elements under diagonal
115 #if 1
116 // evaluate givens rotation
117 valueType theta = evalAngle( n, A, lda, iAbs, jAbs ) ;
118 // apply givens rotation
119 GivensLeft ( n, A, lda, iAbs, jAbs, theta ) ;
120 GivensRight( n, A, lda, iAbs, jAbs, theta ) ;
121 // apply givens rotation to Q
122 GivensRight( n, Q, ldq, iAbs, jAbs, theta ) ;
123 #else
124 for ( indexType i = 1 ; i < n ; ++i ) {
125 for ( indexType j = 0 ; j < i ; ++j ) {
126 // evaluate givens rotation
127 valueType theta = evalAngle( n, A, lda, i, j ) ;
128 // apply givens rotation
129 GivensLeft ( n, A, lda, i, j, theta ) ;
130 GivensRight( n, A, lda, i, j, theta ) ;
131 // apply givens rotation to Q
132 GivensRight( n, Q, ldq, i, j, theta ) ;
133 }
134 }
135 #endif
136 }
137 return false ;
138 }
139}
A driver code
File test_jacobi.cc
1#include "Jacobi.hh"
2#include "alglin.hh"
3#include <iostream>
4#include <iomanip>
5
6using Jacobi::valueType ;
7using Jacobi::indexType ;
8
9using namespace std ;
10
11#define NRow 4
12#define LDA 4
13#define LDQ 4
14#define LDAQ 4
15
16int
17main() {
18
19 valueType Q[NRow*LDQ], AQ[NRow*LDAQ], Asaved[NRow*LDA] ;
20 valueType A[NRow*LDA] = {
21 4, 1, 1, 1,
22 1, 4, 1, 1,
23 1, 1, 4, 1,
24 1, 1, 1, 4
25 } ;
26
27 alglin::copy( NRow*LDA, A, 1, Asaved, 1) ;
28
29 bool ok = Jacobi::Jacobi( NRow, A, LDA, Q, LDQ, 10E-12, 100000 ) ;
30
31 cout << "\nDiagonalized matrix\n\n" ;
32 for ( indexType i = 0 ; i < NRow ; ++i ) {
33 for ( indexType j = 0 ; j < NRow ; ++j ) {
34 cout << setw(10) << A[ i + LDA * j ] << ' ' ;
35 }
36 cout << '\n' ;
37 }
38
39 cout << "\nEigenvectors\n" ;
40 cout << '\n' ;
41
42 for ( indexType i = 0 ; i < NRow ; ++i ) {
43 for ( indexType j = 0 ; j < NRow ; ++j ) {
44 cout << setw(10) << Q[ i + LDQ * j ] << ' ' ;
45 }
46 cout << '\n' ;
47 }
48
49 alglin::gemm( "No Transpose", "No Transpose", NRow, NRow, NRow,
50 1.0, Asaved, LDA,
51 Q, LDQ,
52 0.0, AQ, LDAQ ) ;
53
54 alglin::gemm( "No Transpose", "No Transpose", NRow, NRow, NRow,
55 1.0, Q, LDQ,
56 A, LDA,
57 -1.0, AQ, LDAQ ) ;
58
59 cout << "\nResidual AQ - Q Lambda\n" ;
60 cout << '\n' ;
61
62 for ( indexType i = 0 ; i < NRow ; ++i ) {
63 for ( indexType j = 0 ; j < NRow ; ++j ) {
64 cout << setw(10) << AQ[ i + LDQ * j ] << ' ' ;
65 }
66 cout << '\n' ;
67 }
68
69 return 0 ;
70}