1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #include "lsqrDense.h"
20 #include "vcl_compiler_detection.h"
21
lsqrDense()22 lsqrDense::lsqrDense()
23 {
24 this->A = nullptr;
25 }
26
27
28 lsqrDense::~lsqrDense() = default;
29
30
31 void
SetMatrix(double ** inputA)32 lsqrDense::SetMatrix( double ** inputA )
33 {
34 this->A = inputA;
35 }
36
37
38 /**
39 * computes y = y + A*x without altering x.
40 */
41 void lsqrDense::
Aprod1(unsigned int m,unsigned int n,const double * x,double * y) const42 Aprod1(unsigned int m, unsigned int n, const double * x, double * y ) const
43 {
44 for ( unsigned int row = 0; row < m; row++ )
45 {
46 const double * rowA = this->A[row];
47 double sum = 0.0;
48
49 for ( unsigned int col = 0; col < n; col++ )
50 {
51 sum += rowA[col] * x[col];
52 }
53
54 y[row] += sum;
55 }
56 }
57
58
59 /**
60 * computes x = x + A'*y without altering y.
61 */
62 void lsqrDense::
Aprod2(unsigned int m,unsigned int n,double * x,const double * y) const63 Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const
64 {
65 for ( unsigned int col = 0; col < n; col++ )
66 {
67 double sum = 0.0;
68
69 for ( unsigned int row = 0; row < m; row++ )
70 {
71 sum += this->A[row][col] * y[row];
72 }
73
74 x[col] += sum;
75 }
76 }
77
78
79 /*
80
81 returns x = (I - 2*z*z')*x.
82
83 Implemented as x = x - z * ( 2*( z'*x ) )
84
85 */
86 void lsqrDense::
HouseholderTransformation(unsigned int n,const double * z,double * x) const87 HouseholderTransformation(unsigned int n, const double * z, double * x ) const
88 {
89 // First, compute z'*x as a scalar product.
90 double scalarProduct = 0.0;
91 const double * zp = z;
92 const double * zend = zp + n;
93 double * xp = x;
94 while( zp != zend )
95 {
96 scalarProduct += (*zp++) * (*xp++);
97 }
98
99 // Second, double the value of the scalar product.
100 scalarProduct += scalarProduct;
101
102 // Last, compute x = x - z * (2*z'*x) This subtract from x, double
103 // the componenent of x that is parallel to z, effectively reflecting
104 // x across the hyperplane whose normal is defined by z.
105 zp = z;
106 xp = x;
107 zend = zp + n;
108 while( zp != zend )
109 {
110 *xp++ -= scalarProduct * (*zp++);
111 }
112 }
113