Lesson of 19 February 2010¶
Lesson file in a zip file
The header of QR decomposition
File qr.hh
1#ifndef QR_HH
2#define QR_HH
3
4namespace QR {
5
6 typedef double valueType ; // the type for real number
7 typedef unsigned indexType ; // the type for integer (unsigned)
8
9 indexType
10 factorize( indexType const n, /* Number of rows of matrix A */
11 indexType const m, /* Number of columns of matrix A */
12 valueType A[], /* Pointer to coefficients of A */
13 indexType const lda, /* Leading dimension of A, Fortran style */
14 indexType ipiv[], /* Permutation vector */
15 valueType work[] ) ;
16
17 void
18 solve( indexType const n, /* Number of rows and columns of matrix A */
19 valueType const QR[], /* Pointer to coefficients of LU */
20 indexType const lda, /* Leading dimension of QR, Fortran style */
21 indexType const ipiv[], /* Permutation vector */
22 valueType const work[],
23 valueType b[], /* Pointer to coefficients of vector b */
24 indexType incB ) ;
25
26 void
27 applyReflections( indexType const n, /* Number of rows of matrix QR */
28 indexType const m, /* Number of columns of matrix QR */
29 valueType const QR[], /* Pointer to coefficients of QR */
30 indexType const lda, /* Leading dimension of QR, Fortran style */
31 indexType const ipiv[], /* Permutation vector */
32 valueType const work[],
33 valueType v[],
34 indexType incV ) ;
35
36}
37
38#endif
The code of QR decomposition
File qr.cc
1#include "QR.hh"
2#include "alglin.hh"
3
4namespace QR {
5
6 // Purpose
7 // =======
8 // generates a real elementary reflector H of order n, such that
9 //
10 // H * ( x0 ) = ( beta ), H^T * H = I.
11 // ( v ) ( 0 )
12 //
13 // where alpha and beta are scalars, beta is non-negative, and x is
14 // an (n-1)-element real vector. H is represented in the form
15 //
16 // 2
17 // H = I - --------------- ( alpha ) * ( alpha v' ),
18 // alpha^2 + v^2 ( v )
19 //
20 // 2
21 // H * ( x0 ) = ( x0 ) - -------------- ( alpha ) * ( alpha v' ) ( x0 )
22 // ( v ) ( v ) alpha^2+v^2 ( v ) ( v )
23 //
24 // 2*(x0*alpha+v^2)
25 // = ( x0 ) - ------------------ ( alpha )
26 // ( v ) alpha^2+v^2 ( v )
27 //
28 // alpha = sign(x0) * (|x0| + || (x0,xv) ||)
29 //
30 // 2 alpha * ( alpha x0 + v^2)
31 // beta = x0 - ------------------------------
32 // alpha^2 + v^2
33 //
34 // H = I - gamma ( 1 ) * ( 1 v'/alpha ),
35 // ( v/alpha )
36 //
37 // gamma = 2*alpha^2/(alpha^2+v'*v)
38 //
39
40 valueType // The value tau.
41 buildReflection( indexType n, // The order of the elementary reflector.
42 valueType x[], // On entry, the vector x. On exit, it is overwritten with the vector v and beta
43 valueType xNorm ) { // The norm of vector x
44 // xNorm^2 = x0^2 + v^2
45 valueType x0 = x[0] ;
46 valueType xNorm2 = xNorm*xNorm ;
47 valueType vNorm2 = xNorm2-x0*x0 ;
48 valueType alpha = ( x0 > 0 ? x0 + xNorm : x0 - xNorm ) ;
49 valueType gamma = 2/(1+vNorm2/(alpha*alpha)) ;
50 x[0] -= gamma * ( x0 + vNorm2 / alpha ) ; // beta
51 alglin::scal( n-1, 1/alpha, x+1, 1 ) ;
52 return gamma ;
53 }
54
55 indexType
56 factorize( indexType const n, /* Number of rows of matrix A */
57 indexType const m, /* Number of columns of matrix A */
58 valueType A[], /* Pointer to coefficients of A */
59 indexType const lda, /* Leading dimension of A, Fortran style */
60 indexType ipiv[], /* Permutation vector */
61 valueType work[] ) {
62 indexType nm = m < n ? m : n-1 ;
63 for ( indexType i = 0 ; i < nm ; ++i ) {
64 // evaluate column norm
65 for ( indexType j = i ; j < m ; ++j ) work[j] = alglin::nrm2( n-i, A+i+j*lda, 1 ) ;
66 // find the column with maximum norm
67 indexType ip = i+alglin::iamax( m-i, work+i, 1 ) ;
68 if ( work[ip] == 0 ) return i ; // no more column found
69 // swap columns
70 if ( i != ip ) alglin::swap( n, A + i*lda, 1, A + ip*lda, 1 ) ;
71 ipiv[i] = ip ; // store the pivot index to vector ipiv
72 // build householder reflection
73 valueType * Ai = A + i*(lda+1) ;
74 work[i] = buildReflection( n-i, Ai, work[ip] ) ;
75 // apply householder reflection to submatrix
76 valueType * w = work + i + 1 ;
77 // Ai = ( beta v )^T not (1 v^T)
78 valueType tmp = *Ai ;
79 *Ai = 1 ;
80 // gemv: w <- beta * w + alpha * A^T x
81 // row column alpha A x beta w
82 alglin::gemv( "Transpose", n-i, m-i-1, 1.0, Ai+lda, lda, Ai, 1, 0.0, w, 1 ) ;
83 // ger: A <- A + alpha u v^T
84 // row column alpha A u v
85 alglin::ger( n-i, m-i-1, -work[i], Ai, 1, w, 1, Ai+lda, lda ) ;
86 *Ai = tmp ;
87 }
88 return 0 ;
89 }
90
91 void
92 applyReflections( indexType const n, /* Number of rows of matrix QR */
93 indexType const m, /* Number of columns of matrix QR */
94 valueType const QR[], /* Pointer to coefficients of QR */
95 indexType const lda, /* Leading dimension of QR, Fortran style */
96 indexType const ipiv[], /* Permutation vector */
97 valueType const work[],
98 valueType v[],
99 indexType incV ) {
100 indexType nm = m < n ? m : n-1 ;
101 for ( indexType i = 0 ; i < nm ; ++i ) {
102 valueType * vi = v + i*incV ;
103 valueType const * Qi = QR + i*lda + i + 1 ;
104 valueType bf = work[i] * (*vi + alglin::dot( n-i-1, Qi, 1, vi + incV, incV )) ;
105 *vi -= bf ;
106 alglin::axpy( n-i-1, -bf, Qi, 1, vi + incV, 1 ) ;
107 }
108 }
109
110 void
111 solve( indexType const n, /* Number of rows and columns of matrix A */
112 valueType const QR[], /* Pointer to coefficients of LU */
113 indexType const lda, /* Leading dimension of LU, Fortran style */
114 indexType const ipiv[], /* Permutation vector */
115 valueType const work[],
116 valueType b[], /* Pointer to coefficients of vector b */
117 indexType incB ) {
118 applyReflections( n, n, QR, lda, ipiv, work, b, incB ) ;
119 alglin::trsv( "Upper", "No Transpose", "Non Unit", n, QR, lda, b, incB ) ;
120 for ( indexType i = 0 ; i < n-1 ; ++i ) {
121 indexType ip = ipiv[i] ;
122 if ( ip != i ) std::swap( b[i*incB], b[ip*incB] ) ;
123 }
124
125 }
126}