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 "linalg/lsqrDense.h"
20
21 #include <iostream>
22
23 #include <cstdlib>
24 #include <cmath>
25
26
main(int,char * [])27 int main( int , char * [] )
28 {
29
30 lsqrDense solver;
31
32 const double tolerance = 1e-9;
33
34 const unsigned int n = 2;
35 double x[n];
36 double z[n];
37
38 x[0] = 3.0;
39 x[1] = 5.0;
40
41 z[0] = 0.0;
42 z[1] = 1.0;
43
44 solver.HouseholderTransformation(n,z,x);
45
46 std::cout << x[0] << " " << x[1] << std::endl;
47
48 { // Test 1 Dnrm2()
49 const unsigned int n1 = 5;
50 double x1[n1];
51 x1[0] = 1.0;
52 x1[1] = 1.0;
53 x1[2] = 1.0;
54 x1[3] = 1.0;
55 x1[4] = 1.0;
56
57 const double norm =solver.Dnrm2( n1, x1 );
58 const double expectedNorm = sqrt(5.0);
59
60 const double ratioOfDifference =
61 fabs( norm - expectedNorm ) / expectedNorm;
62
63 if( ratioOfDifference > tolerance )
64 {
65 std::cerr << "Error in Dnrm2() test 1 " << std::endl;
66 std::cerr << "Expected = " << expectedNorm << std::endl;
67 std::cerr << "Received = " << norm << std::endl;
68 std::cerr << "ratioOfDifference = " << ratioOfDifference << std::endl;
69 return EXIT_FAILURE;
70 }
71 std::cout << "Dnrm2 test 1 passed " << std::endl;
72 }
73
74 { // Test 2 Dnrm2()
75 const unsigned int n2 = 5;
76
77 const double dominantValue = 1e+300;
78
79 double x2[n2];
80 x2[0] = 1e+30;
81 x2[1] = 1e+200;
82 x2[2] = dominantValue;
83 x2[3] = 1e+2;
84 x2[4] = 1e+1;
85
86 const double norm =solver.Dnrm2( n2, x2 );
87 const double expectedNorm = dominantValue;
88
89 const double ratioOfDifference =
90 fabs( norm - expectedNorm ) / expectedNorm;
91
92 if( ratioOfDifference > tolerance )
93 {
94 std::cerr << "Error in Dnrm2() test 2 " << std::endl;
95 std::cerr << "Expected = " << expectedNorm << std::endl;
96 std::cerr << "Received = " << norm << std::endl;
97 std::cerr << "ratioOfDifference = " << ratioOfDifference << std::endl;
98 return EXIT_FAILURE;
99 }
100 std::cout << "Dnrm2 test 2 passed " << std::endl;
101 }
102
103
104 { // Testing D2Norm
105 const double dominantValue = 1e+300;
106 const double minorValue = 1e-100;
107 const double a = dominantValue;
108 const double b = minorValue;
109 const double d2norm = solver.D2Norm( a, b );
110
111 const double ratioOfDifference =
112 fabs( d2norm - dominantValue ) / dominantValue;
113
114 if( ratioOfDifference > tolerance )
115 {
116 std::cerr << "Error in D2Norm() test 1 " << std::endl;
117 std::cerr << "Expected = " << dominantValue << std::endl;
118 std::cerr << "Received = " << d2norm << std::endl;
119 return EXIT_FAILURE;
120 }
121 std::cout << "D2Norm test 1 passed " << std::endl;
122 }
123
124 { // Testing D2Norm
125 const double friendlyValue1 = 1e+3;
126 const double friendlyValue2 = 1e+2;
127 const double a = friendlyValue1;
128 const double b = friendlyValue2;
129 const double d2norm = solver.D2Norm( a, b );
130 const double expectedD2norm = sqrt( a*a + b*b );
131
132 const double ratioOfDifference =
133 fabs( d2norm - expectedD2norm ) / expectedD2norm;
134
135 if( ratioOfDifference > tolerance )
136 {
137 std::cerr << "Error in D2Norm() test 2 " << std::endl;
138 std::cerr << "Expected = " << expectedD2norm << std::endl;
139 std::cerr << "Received = " << d2norm << std::endl;
140 return EXIT_FAILURE;
141 }
142 std::cout << "D2Norm test 2 passed " << std::endl;
143 }
144
145
146 { // Testing D2Norm
147 const double zero = 0.0;
148 const double a = zero;
149 const double b = zero;
150 const double d2norm = solver.D2Norm( a, b );
151
152 const double difference = fabs( d2norm - zero );
153
154 if( difference > tolerance )
155 {
156 std::cerr << "Error in D2Norm() test 3 " << std::endl;
157 std::cerr << "Expected = " << zero << std::endl;
158 std::cerr << "Received = " << d2norm << std::endl;
159 return EXIT_FAILURE;
160 }
161 std::cout << "D2Norm test 3 passed " << std::endl;
162 }
163
164
165 solver.SetOutputStream( std::cout );
166
167 const double eps = 1e-15;
168
169 solver.SetEpsilon( eps );
170
171 const double damp = 0.0;
172
173 solver.SetDamp( damp );
174
175 solver.SetMaximumNumberOfIterations( 100 );
176
177 solver.SetToleranceA( 1e-16 );
178 solver.SetToleranceB( 1e-16 );
179
180 solver.SetUpperLimitOnConditional( 1.0 / ( 10 * sqrt( eps ) ) );
181
182 solver.SetStandardErrorEstimatesFlag( true );
183
184
185 { // basic test for Scale()
186 const double factor = 8.0;
187 const unsigned int mm = 5;
188 double xx[mm];
189 for ( unsigned int i = 0; i < mm; i++ )
190 {
191 xx[i] = i * 5.0;
192 }
193
194 solver.Scale( mm, factor, xx );
195
196 for ( unsigned int i = 0; i < mm; i++ )
197 {
198 const double expectedValue = ( i*5.0*factor );
199 const double ratioOfDifference =
200 fabs( xx[i] - expectedValue ) / expectedValue;
201
202 if ( ratioOfDifference > tolerance )
203 {
204 std::cerr << "Error in method Scale() " << std::endl;
205 return EXIT_FAILURE;
206 }
207 }
208 }
209
210
211 { // basic test for Aprod1()
212 const unsigned int mm = 2;
213 const unsigned int nn = 2;
214 double xx[nn];
215 xx[0] = 1.0;
216 xx[1] = 7.0;
217 using RowType = double *;
218 RowType A[mm];
219 double AA[4];
220 A[0] = &(AA[0]);
221 A[1] = &(AA[2]);
222 A[0][0] = 1.0;
223 A[0][1] = 0.0;
224 A[1][0] = 0.0;
225 A[1][1] = 1.0;
226 solver.SetMatrix( A );
227 double yy[mm];
228 yy[0] = 0.0;
229 yy[1] = 0.0;
230 solver.Aprod1( mm, nn, xx, yy );
231 std::cout << "yy = " << yy[0] << " " << yy[1] << std::endl;
232 }
233
234 { // basic test for Aprod1()
235 const unsigned int mm = 2;
236 const unsigned int nn = 3;
237 double xx[nn];
238 xx[0] = 1.0;
239 xx[1] = 7.0;
240 xx[2] = 9.0;
241 using RowType = double *;
242 RowType A[mm];
243 double AA[6];
244 A[0] = &(AA[0]);
245 A[1] = &(AA[3]);
246 A[0][0] = 1.0;
247 A[0][1] = 2.0;
248 A[0][2] = 3.0;
249 A[1][0] = 4.0;
250 A[1][1] = 5.0;
251 A[1][2] = 6.0;
252 solver.SetMatrix( A );
253 double yy[mm];
254 yy[0] = 0.0;
255 yy[1] = 0.0;
256 solver.Aprod1( mm, nn, xx, yy );
257 std::cout << "yy = " << yy[0] << " " << yy[1] << std::endl;
258 }
259
260 { // basic test for Aprod1()
261 const unsigned int mm = 3;
262 const unsigned int nn = 2;
263 double xx[nn];
264 xx[0] = 1.0;
265 xx[1] = 7.0;
266 using RowType = double *;
267 RowType A[mm];
268 double AA[6];
269 A[0] = &(AA[0]);
270 A[1] = &(AA[2]);
271 A[2] = &(AA[4]);
272 A[0][0] = 1.0;
273 A[0][1] = 2.0;
274 A[1][0] = 3.0;
275 A[1][1] = 4.0;
276 A[2][0] = 5.0;
277 A[2][1] = 6.0;
278 solver.SetMatrix( A );
279 double yy[mm];
280 yy[0] = 0.0;
281 yy[1] = 0.0;
282 yy[2] = 0.0;
283 solver.Aprod1( mm, nn, xx, yy );
284 std::cout << "yy = " << yy[0] << " " << yy[1] << " " << yy[2] << std::endl;
285 }
286
287 { // basic test for Aprod2()
288 const unsigned int mm = 2;
289 const unsigned int nn = 2;
290 double xx[nn];
291 xx[0] = 0.0;
292 xx[1] = 0.0;
293 using RowType = double *;
294 RowType A[mm];
295 double AA[4];
296 A[0] = &(AA[0]);
297 A[1] = &(AA[2]);
298 A[0][0] = 1.0;
299 A[0][1] = 0.0;
300 A[1][0] = 0.0;
301 A[1][1] = 1.0;
302 solver.SetMatrix( A );
303 double yy[mm];
304 yy[0] = 1.0;
305 yy[1] = 7.0;
306 solver.Aprod2( mm, nn, xx, yy );
307 std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl;
308 }
309
310 { // basic test for Aprod2()
311 const unsigned int mm = 2;
312 const unsigned int nn = 3;
313 double xx[nn];
314 xx[0] = 0.0;
315 xx[1] = 0.0;
316 xx[2] = 0.0;
317 using RowType = double *;
318 RowType A[mm];
319 double AA[6];
320 A[0] = &(AA[0]);
321 A[1] = &(AA[3]);
322 A[0][0] = 1.0;
323 A[0][1] = 2.0;
324 A[0][2] = 3.0;
325 A[1][0] = 4.0;
326 A[1][1] = 5.0;
327 A[1][2] = 6.0;
328 solver.SetMatrix( A );
329 double yy[mm];
330 yy[0] = 1.0;
331 yy[1] = 5.0;
332 solver.Aprod2( mm, nn, xx, yy );
333 std::cout << "xx = " << xx[0] << " " << xx[1] << " " << xx[2] << std::endl;
334 }
335
336 { // basic test for Aprod2()
337 const unsigned int mm = 3;
338 const unsigned int nn = 2;
339 double xx[nn];
340 xx[0] = 0.0;
341 xx[1] = 0.0;
342 using RowType = double *;
343 RowType A[mm];
344 double AA[6];
345 A[0] = &(AA[0]);
346 A[1] = &(AA[2]);
347 A[2] = &(AA[4]);
348 A[0][0] = 1.0;
349 A[0][1] = 2.0;
350 A[1][0] = 3.0;
351 A[1][1] = 4.0;
352 A[2][0] = 5.0;
353 A[2][1] = 6.0;
354 solver.SetMatrix( A );
355 double yy[mm];
356 yy[0] = 1.0;
357 yy[1] = 7.0;
358 yy[2] = 9.0;
359 solver.Aprod2( mm, nn, xx, yy );
360 std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl;
361 }
362
363 { // basic test for Solve()
364 const unsigned int mm = 2;
365 const unsigned int nn = 2;
366 double bb[nn];
367 bb[0] = 1.0;
368 bb[1] = 7.0;
369 double xx[mm];
370 solver.SetStandardErrorEstimatesFlag( false );
371 using RowType = double *;
372 RowType A[mm];
373 double AA[4];
374 A[0] = &(AA[0]);
375 A[1] = &(AA[2]);
376 A[0][0] = 1.0;
377 A[0][1] = 0.0;
378 A[1][0] = 0.0;
379 A[1][1] = 1.0;
380 solver.SetMatrix( A );
381 solver.SetMaximumNumberOfIterations( 10 );
382 solver.Solve( mm, nn, bb, xx );
383 std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl;
384 }
385
386 { // basic test for Solve()
387 const unsigned int mm = 2;
388 const unsigned int nn = 2;
389 double bb[nn];
390 bb[0] = 1.0;
391 bb[1] = 7.0;
392 double xx[mm];
393 solver.SetStandardErrorEstimatesFlag( true );
394 double se[nn];
395 solver.SetStandardErrorEstimates( se );
396 using RowType = double *;
397 RowType A[mm];
398 double AA[4];
399 A[0] = &(AA[0]);
400 A[1] = &(AA[2]);
401 A[0][0] = 1.0;
402 A[0][1] = 0.0;
403 A[1][0] = 0.0;
404 A[1][1] = 1.0;
405 solver.SetMatrix( A );
406 solver.SetMaximumNumberOfIterations( 10 );
407 solver.Solve( mm, nn, bb, xx );
408 std::cout << "se = " << se[0] << " " << se[1] << std::endl;
409 std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl;
410 }
411
412
413 { // basic test for Solve()
414 const unsigned int mm = 2;
415 const unsigned int nn = 2;
416 double bb[nn];
417 bb[0] = -9.0;
418 bb[1] = 15.0;
419 double xx[mm];
420 solver.SetStandardErrorEstimatesFlag( true );
421 double se[nn];
422 solver.SetStandardErrorEstimates( se );
423 using RowType = double *;
424 RowType A[mm];
425 double AA[4];
426 A[0] = &(AA[0]);
427 A[1] = &(AA[2]);
428 A[0][0] = 1.0;
429 A[0][1] = 0.0;
430 A[1][0] = 0.0;
431 A[1][1] = 1.0;
432 solver.SetMatrix( A );
433 solver.SetMaximumNumberOfIterations( 10 );
434 solver.Solve( mm, nn, bb, xx );
435 std::cout << "se = " << se[0] << " " << se[1] << std::endl;
436 std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl;
437 }
438
439
440 std::cout << "Stopped because " << solver.GetStoppingReason() << std::endl;
441 std::cout << "Used " << solver.GetNumberOfIterationsPerformed() << " Iterations" << std::endl;
442 std::cout << "Frobenius norm estimation of Abar = " << solver.GetFrobeniusNormEstimateOfAbar() << std::endl;
443 std::cout << "Condition number estimation of Abar = " << solver.GetConditionNumberEstimateOfAbar() << std::endl;
444 std::cout << "Estimate of final value of norm(rbar) = " << solver.GetFinalEstimateOfNormRbar() << std::endl;
445 std::cout << "Estimate of final value of norm of residuals = " << solver.GetFinalEstimateOfNormOfResiduals() << std::endl;
446 std::cout << "Estimate of norm of final solution = " << solver.GetFinalEstimateOfNormOfX() << std::endl;
447
448 return EXIT_SUCCESS;
449 }
450