1 # include <cstdlib>
2 # include <iostream>
3 # include <iomanip>
4 # include <cstring>
5
6 # include "sandia_rules.hpp"
7 # include "sandia_sgmgg.hpp"
8
9 int main ( );
10 void test01 ( );
11 void test02 ( );
12 void test03 ( );
13 void test04 ( );
14 void test05 ( );
15 void test06 ( );
16 void sandia_sgmgg_coef_naive_test ( int dim_num, int point_num,
17 int sparse_index[] );
18 void sandia_sgmgg_neighbor_naive_test ( int dim_num, int point_num,
19 int sparse_index[] );
20
21 //****************************************************************************80
22
main()23 int main ( )
24
25 //****************************************************************************80
26 //
27 // Purpose:
28 //
29 // SANDIA_SGMGG_PRB tests generalized sparse grid routines.
30 //
31 // Licensing:
32 //
33 // This code is distributed under the GNU LGPL license.
34 //
35 // Modified:
36 //
37 // 23 August 2011
38 //
39 // Author:
40 //
41 // John Burkardt
42 //
43 {
44 webbur::timestamp ( );
45 std::cout << "\n";
46 std::cout << "SANDIA_SGMGG_PRB:\n";
47 std::cout << " C++ version.\n";
48 std::cout << " Test the SANDIA_SGMGG library.\n";
49
50 test01 ( );
51 test02 ( );
52 test03 ( );
53 test04 ( );
54 test05 ( );
55 test06 ( );
56 //
57 // Terminate.
58 //
59 std::cout << "\n";
60 std::cout << "SANDIA_SGMGG_PRB:\n";
61 std::cout << " Normal end of execution.\n";
62 std::cout << "\n";
63 webbur::timestamp ( );
64
65 return 0;
66 }
67 //****************************************************************************80
68
test01()69 void test01 ( )
70
71 //****************************************************************************80
72 //
73 // Purpose:
74 //
75 // TEST01 demonstrates naive coefficient calculation.
76 //
77 // Licensing:
78 //
79 // This code is distributed under the GNU LGPL license.
80 //
81 // Modified:
82 //
83 // 05 August 2010
84 //
85 // Author:
86 //
87 // John Burkardt
88 //
89 {
90 int dim_num;
91 int point_num;
92 int sparse_index1[2*7] = {
93 0, 2,
94 0, 3,
95 1, 1,
96 1, 2,
97 2, 0,
98 2, 1,
99 3, 0 };
100 int sparse_index2[3*19] = {
101 0, 1, 0,
102 0, 2, 0,
103 0, 3, 0,
104 1, 0, 0,
105 1, 1, 0,
106 1, 2, 0,
107 2, 0, 0,
108 2, 1, 0,
109 3, 0, 0,
110 0, 0, 1,
111 0, 1, 1,
112 0, 2, 1,
113 1, 0, 1,
114 1, 1, 1,
115 2, 0, 1,
116 0, 0, 2,
117 0, 1, 2,
118 1, 0, 2,
119 0, 0, 3 };
120 int sparse_index3[2*8] = {
121 0, 2,
122 1, 1,
123 1, 2,
124 2, 1,
125 3, 0,
126 3, 1,
127 4, 0,
128 5, 0 };
129 int sparse_index4[2*8] = {
130 0, 0,
131 0, 1,
132 0, 2,
133 0, 3,
134 1, 0,
135 1, 1,
136 2, 0,
137 3, 0 };
138
139 std::cout << "\n";
140 std::cout << "TEST01:\n";
141 std::cout << " Demonstrate naive coefficient calculations.\n";
142 //
143 // Isotropic grid in 2D.
144 //
145 std::cout << "\n";
146 std::cout << " 1) Isotropic grid in 2D:\n";
147 dim_num = 2;
148 point_num = 7;
149
150 sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index1 );
151 //
152 // Isotropic grid in 3D.
153 //
154 std::cout << "\n";
155 std::cout << " 2) Isotropic grid in 3D:\n";
156 dim_num = 3;
157 point_num = 19;
158
159 sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index2 );
160 //
161 // Anisotropic grid in 2D.
162 //
163 std::cout << "\n";
164 std::cout << " 3) Ansotropic grid in 2D:\n";
165 dim_num = 2;
166 point_num = 8;
167
168 sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index3 );
169 //
170 // Generalized grid in 2D.
171 //
172 std::cout << "\n";
173 std::cout << " 4) Generalized grid in 2D:\n";
174 dim_num = 2;
175 point_num = 8;
176
177 sandia_sgmgg_coef_naive_test ( dim_num, point_num, sparse_index4 );
178
179 return;
180 }
181 //****************************************************************************80
182
test02()183 void test02 ( )
184
185 //****************************************************************************80
186 //
187 // Purpose:
188 //
189 // TEST02 demonstrates naive neighbor calculation.
190 //
191 // Licensing:
192 //
193 // This code is distributed under the GNU LGPL license.
194 //
195 // Modified:
196 //
197 // 06 August 2010
198 //
199 // Author:
200 //
201 // John Burkardt
202 //
203 {
204 int dim_num;
205 int point_num;
206 int *sparse_index;
207 int sparse_index1[2*7] = {
208 0, 2,
209 0, 3,
210 1, 1,
211 1, 2,
212 2, 0,
213 2, 1,
214 3, 0 };
215 int sparse_index2[3*19] = {
216 0, 1, 0,
217 0, 2, 0,
218 0, 3, 0,
219 1, 0, 0,
220 1, 1, 0,
221 1, 2, 0,
222 2, 0, 0,
223 2, 1, 0,
224 3, 0, 0,
225 0, 0, 1,
226 0, 1, 1,
227 0, 2, 1,
228 1, 0, 1,
229 1, 1, 1,
230 2, 0, 1,
231 0, 0, 2,
232 0, 1, 2,
233 1, 0, 2,
234 0, 0, 3 };
235 int sparse_index3[2*8] = {
236 0, 2,
237 1, 1,
238 1, 2,
239 2, 1,
240 3, 0,
241 3, 1,
242 4, 0,
243 5, 0 };
244 int sparse_index4[2*8] = {
245 0, 0,
246 0, 1,
247 0, 2,
248 0, 3,
249 1, 0,
250 1, 1,
251 2, 0,
252 3, 0 };
253
254 std::cout << "\n";
255 std::cout << "TEST02:\n";
256 std::cout << " Demonstrate naive neighbor calculations.\n";
257 //
258 // Isotropic grid in 2D.
259 //
260 std::cout << "\n";
261 std::cout << " 1) Isotropic grid in 2D:\n";
262 dim_num = 2;
263 point_num = 7;
264
265 sandia_sgmgg_neighbor_naive_test ( dim_num, point_num, sparse_index1 );
266 //
267 // Isotropic grid in 3D.
268 //
269 std::cout << "\n";
270 std::cout << " 2) Isotropic grid in 3D:\n";
271 dim_num = 3;
272 point_num = 19;
273
274 sandia_sgmgg_neighbor_naive_test ( dim_num, point_num, sparse_index2 );
275 //
276 // Anisotropic grid in 2D.
277 //
278 std::cout << "\n";
279 std::cout << " 3) Ansotropic grid in 2D:\n";
280 dim_num = 2;
281 point_num = 8;
282
283 sandia_sgmgg_neighbor_naive_test ( dim_num, point_num, sparse_index3 );
284 //
285 // Generalized grid in 2D.
286 //
287 std::cout << "\n";
288 std::cout << " 4) Generalized grid in 2D:\n";
289 dim_num = 2;
290 point_num = 8;
291
292 sandia_sgmgg_neighbor_naive_test ( dim_num, point_num, sparse_index4 );
293
294 return;
295 }
296 //****************************************************************************80
297
test03()298 void test03 ( )
299
300 //****************************************************************************80
301 //
302 // Purpose:
303 //
304 // TEST03 sets up the GG data structure.
305 //
306 // Licensing:
307 //
308 // This code is distributed under the GNU LGPL license.
309 //
310 // Modified:
311 //
312 // 20 August 2010
313 //
314 // Author:
315 //
316 // John Burkardt
317 //
318 {
319 # define GG_NA_01 4
320 # define GG_ND_01 2
321 # define GG_NI_01 10
322 # define GG_NO_01 6
323
324 int *gg_a;
325 int gg_a_01[GG_NA_01] = { 3, 6, 8, 9 };
326 int *gg_b;
327 int gg_b_01[GG_ND_01*GG_NI_01] = {
328 -1, -1,
329 -1, 0,
330 -1, 1,
331 -1, 2,
332 0, -1,
333 1, 4,
334 2, 5,
335 4, -1,
336 5, 7,
337 7, -1 };
338 int *gg_f;
339 int gg_f_01[GG_ND_01*GG_NI_01] = {
340 4, 1,
341 5, 2,
342 6, 3,
343 -1, -1,
344 7, 5,
345 8, 6,
346 -1, -1,
347 9, 8,
348 -1, -1,
349 -1, -1 };
350 double *gg_g;
351 double gg_g_01[GG_NI_01] = {
352 0.1, 1.1, 2.2, 3.0, 1.0,
353 2.1, 3.2, 2.0, 3.3, 3.1 };
354 int *gg_i;
355 int gg_i_01[GG_ND_01 * GG_NI_01] = {
356 0, 0,
357 0, 1,
358 0, 2,
359 0, 3,
360 1, 0,
361 1, 1,
362 1, 2,
363 2, 0,
364 2, 1,
365 3, 0 };
366 int gg_ma = 100;
367 int gg_mi = 100;
368 int gg_mo = 100;
369 int gg_na;
370 int gg_nd;
371 int gg_ni;
372 int gg_no;
373 int *gg_o;
374 int gg_o_01[GG_NO_01] = { 0, 1, 2, 4, 5, 7 };
375 int i;
376 int max_index;
377
378 std::cout << "\n";
379 std::cout << "TEST03:\n";
380 std::cout << " Set up examples of the GG (Gerstner-Griebel)\n";
381 std::cout << " data structure for generalized sparse grids.\n";
382 //
383 // Isotropic grid in 2D.
384 //
385 // 3 | 4
386 // 2 | 3 7
387 // 1 | 2 6 9
388 // 0 | 1 5 8 10
389 // +-----------
390 // 0 1 2 3
391 //
392 std::cout << "\n";
393 std::cout << " 1) Isotropic grid in 2D\n";
394 //
395 // Get sizes for this problem.
396 //
397 gg_na = GG_NA_01;
398 gg_nd = GG_ND_01;
399 gg_ni = GG_NI_01;
400 gg_no = GG_NO_01;
401 //
402 // Create work arrays.
403 // Because these will grow and shrink during the computation, we must give
404 // them an initial size big enough that we don't need to resize.
405 //
406 gg_a = new int [ gg_ma ];
407 gg_b = new int [ gg_nd * gg_mi ];
408 gg_f = new int [ gg_nd * gg_mi ];
409 gg_g = new double [ gg_mi ];
410 gg_i = new int [ gg_nd * gg_mi ];
411 gg_o = new int [ gg_mo ];
412 //
413 // Copy the initial problem data vectors.
414 //
415 webbur::i4vec_copy ( gg_na, gg_a_01, gg_a );
416 webbur::i4mat_copy ( gg_nd, gg_ni, gg_b_01, gg_b );
417 webbur::i4mat_copy ( gg_nd, gg_ni, gg_f_01, gg_f );
418 webbur::r8vec_copy ( gg_ni, gg_g_01, gg_g );
419 webbur::i4mat_copy ( gg_nd, gg_ni, gg_i_01, gg_i );
420 webbur::i4vec_copy ( gg_no, gg_o_01, gg_o );
421 //
422 // Implicit decreasing heap sort on GG_G and GG_A.
423 //
424 std::cout << "\n";
425 std::cout << " Before Heap:\n";
426 std::cout << " I A G\n";
427 std::cout << "\n";
428
429 for ( i = 0; i < gg_na; i++ )
430 {
431 std::cout << " " << std::setw(4) << i
432 << " " << std::setw(4) << gg_a[i]
433 << " " << std::setw(14) << gg_g[gg_a[i]] << "\n";
434 }
435
436 webbur::r8vec_indexed_heap_d ( gg_na, gg_g, gg_a );
437
438 std::cout << "\n";
439 std::cout << " After Heap:\n";
440 std::cout << " I A G\n";
441 std::cout << "\n";
442
443 for ( i = 0; i < gg_na; i++ )
444 {
445 std::cout << " " << std::setw(4) << i
446 << " " << std::setw(4) << gg_a[i]
447 << " " << std::setw(14) << gg_g[gg_a[i]] << "\n";
448 }
449 //
450 // One step of the algorithm is to identify MAX_INDEX, the active index with
451 // maximum G value, and move it to the OLD array.
452 //
453 max_index = webbur::r8vec_indexed_heap_d_extract ( &gg_na, gg_g, gg_a );
454 std::cout << "\n";
455 std::cout << " Transferring index " << max_index
456 << " from active to old set.\n";
457 //
458 // GG_NA was decremented.
459 // The heap structure of GG_A was restored.
460 // To complete the update, we must
461 // append MAX_INDEX to GG_O.
462 // update the active set with the forward neighbors of MAX_INDEX.
463 //
464 gg_o[gg_no] = max_index;
465 gg_no = gg_no + 1;
466 //
467 // Compute forward neighbors of MAX_INDEX.
468 //
469
470 //
471 // DISPLAY CURRENT DATA.
472 //
473 delete [] gg_a;
474 delete [] gg_b;
475 delete [] gg_f;
476 delete [] gg_g;
477 delete [] gg_i;
478 delete [] gg_o;
479
480 return;
481 # undef GG_NA_01
482 # undef GG_ND_01
483 # undef GG_NI_01
484 # undef GG_NO_01
485 }
486 //****************************************************************************80
487
test04()488 void test04 ( )
489
490 //****************************************************************************80
491 //
492 // Purpose:
493 //
494 // TEST04 simulates a set of incremental coefficient calculations.
495 //
496 // Licensing:
497 //
498 // This code is distributed under the GNU LGPL license.
499 //
500 // Modified:
501 //
502 // 11 January 2011
503 //
504 // Author:
505 //
506 // John Burkardt
507 //
508 {
509 int dim_num;
510 int point_num;
511 int point_num2;
512 int sparse_index[2*8] = {
513 0, 0,
514 0, 1,
515 0, 2,
516 0, 3,
517 1, 0,
518 1, 1,
519 2, 0,
520 3, 0 };
521
522 std::cout << "\n";
523 std::cout << "TEST04:\n";
524 std::cout << " Simulate incremental coefficient calculations.\n";
525 //
526 // Generalized grid in 2D.
527 //
528 std::cout << "\n";
529 std::cout << " Generalized grid in 2D:\n";
530 dim_num = 2;
531 point_num = 8;
532
533 for ( point_num2 = 1; point_num2 <= point_num; point_num2++ )
534 {
535 sandia_sgmgg_coef_naive_test ( dim_num, point_num2, sparse_index );
536 }
537 return;
538 }
539 //****************************************************************************80
540
test05()541 void test05 ( )
542
543 //****************************************************************************80
544 //
545 // Purpose:
546 //
547 // TEST05 demonstrates SANDIA_SGMGG_COEF_INC2.
548 //
549 // Licensing:
550 //
551 // This code is distributed under the GNU LGPL license.
552 //
553 // Modified:
554 //
555 // 14 January 2011
556 //
557 // Author:
558 //
559 // John Burkardt
560 //
561 {
562 # define M 2
563 # define N1 5
564
565 int c1[N1] = {
566 +1,
567 +1,
568 +1,
569 -1,
570 -1
571 };
572 int *c3;
573 int i;
574 int j;
575 int m = M;
576 int n1 = N1;
577 int s1[M*N1] = {
578 0, 2,
579 1, 1,
580 2, 0,
581 0, 1,
582 1, 0 };
583 int s2[M] = {
584 1, 2 };
585
586 std::cout << "\n";
587 std::cout << "TEST05:\n";
588 std::cout << " Predict new coefficients given candidate index.\n";
589
590 std::cout << "\n";
591 std::cout << " Spatial dimension M = " << m << "\n";
592 std::cout << " Number of items in active set N1 = " << n1 << "\n";
593 std::cout << "\n";
594 std::cout << " Index Coef Indices\n";
595 std::cout << "\n";
596 for ( j = 0; j < n1; j++ )
597 {
598 std::cout << " " << std::setw(2) << j << ":"
599 << " " << std::setw(4) << c1[j] << " ";
600 for ( i = 0; i < m; i++ )
601 {
602 std::cout << " " << std::setw(2) << s1[i+j*m];
603 }
604 std::cout << "\n";
605 }
606 std::cout << "\n";
607 std::cout << " Candidate: ";
608 for ( i = 0; i < m; i++ )
609 {
610 std::cout << " " << std::setw(2) << s2[i];
611 }
612 std::cout << "\n";
613 //
614 // Generalized grid in 2D.
615 //
616 c3 = new int[n1+1];
617
618 webbur::sandia_sgmgg_coef_inc2 ( m, n1, s1, c1, s2, c3 );
619
620 std::cout << "\n";
621 std::cout << " Index Coef Coef\n";
622 std::cout << " Old New\n";
623 std::cout << "\n";
624 for ( i = 0; i < n1; i++ )
625 {
626 std::cout << " " << std::setw(2) << i << ":"
627 << " " << std::setw(4) << c1[i]
628 << " " << std::setw(4) << c3[i] << "\n";
629 }
630
631 std::cout << " " << std::setw(2) << n1 << ":"
632 << " " << " "
633 << " " << std::setw(4) << c3[n1] << "\n";
634
635 std::cout << " -- ---- ----\n";
636 std::cout << " Sum:"
637 << " " << std::setw(4) << webbur::i4vec_sum ( n1, c1 )
638 << " " << std::setw(4) << webbur::i4vec_sum ( n1 + 1, c3 ) << "\n";
639 delete [] c3;
640
641 return;
642 # undef M
643 # undef N1
644 }
645 //****************************************************************************80
646
test06()647 void test06 ( )
648
649 //****************************************************************************80
650 //
651 // Purpose:
652 //
653 // TEST06 tests a special case for SANDIA_SGMGG_COEF_INC2.
654 //
655 // Licensing:
656 //
657 // This code is distributed under the GNU LGPL license.
658 //
659 // Modified:
660 //
661 // 05 September 2011
662 //
663 // Author:
664 //
665 // John Burkardt
666 //
667 {
668 # define M 2
669 # define N1 5
670
671 int c1[N1] = {
672 +1,
673 +1,
674 +1,
675 -1,
676 -1
677 };
678 int *c3;
679 int i;
680 int j;
681 int m = M;
682 int n1 = N1;
683 int s1[M*N1] = {
684 2, 0,
685 1, 1,
686 0, 2,
687 1, 0,
688 0, 1 };
689 int s2[M] = {
690 3, 0 };
691
692 std::cout << "\n";
693 std::cout << "TEST06:\n";
694 std::cout << " Predict new coefficients given candidate index.\n";
695
696 std::cout << "\n";
697 std::cout << " Spatial dimension M = " << m << "\n";
698 std::cout << " Number of items in active set N1 = " << n1 << "\n";
699 std::cout << "\n";
700 std::cout << " Index Coef Indices\n";
701 std::cout << "\n";
702 for ( j = 0; j < n1; j++ )
703 {
704 std::cout << " " << std::setw(2) << j << ":"
705 << " " << std::setw(4) << c1[j] << " ";
706 for ( i = 0; i < m; i++ )
707 {
708 std::cout << " " << std::setw(2) << s1[i+j*m];
709 }
710 std::cout << "\n";
711 }
712 std::cout << "\n";
713 std::cout << " Candidate: ";
714 for ( i = 0; i < m; i++ )
715 {
716 std::cout << " " << std::setw(2) << s2[i];
717 }
718 std::cout << "\n";
719 //
720 // Generalized grid in 2D.
721 //
722 c3 = new int[n1+1];
723
724 webbur::sandia_sgmgg_coef_inc2 ( m, n1, s1, c1, s2, c3 );
725
726 std::cout << "\n";
727 std::cout << " Index Coef Coef\n";
728 std::cout << " Old New\n";
729 std::cout << "\n";
730 for ( i = 0; i < n1; i++ )
731 {
732 std::cout << " " << std::setw(2) << i << ":"
733 << " " << std::setw(4) << c1[i]
734 << " " << std::setw(4) << c3[i] << "\n";
735 }
736
737 std::cout << " " << std::setw(2) << n1 << ":"
738 << " " << " "
739 << " " << std::setw(4) << c3[n1] << "\n";
740
741 std::cout << " -- ---- ----\n";
742 std::cout << " Sum:"
743 << " " << std::setw(4) << webbur::i4vec_sum ( n1, c1 )
744 << " " << std::setw(4) << webbur::i4vec_sum ( n1 + 1, c3 ) << "\n";
745 delete [] c3;
746
747 return;
748 # undef M
749 # undef N1
750 }
751 //****************************************************************************80
752
sandia_sgmgg_coef_naive_test(int dim_num,int point_num,int sparse_index[])753 void sandia_sgmgg_coef_naive_test ( int dim_num, int point_num,
754 int sparse_index[] )
755
756 //****************************************************************************80
757 //
758 // Purpose:
759 //
760 // SANDIA_SGMGG_COEF_NAIVE_TEST demonstrates "naive" coefficient computations.
761 //
762 // Licensing:
763 //
764 // This code is distributed under the GNU LGPL license.
765 //
766 // Modified:
767 //
768 // 05 August 2010
769 //
770 // Author:
771 //
772 // John Burkardt
773 //
774 {
775 int *coef;
776 int coef_sum;
777
778 webbur::i4mat_transpose_print ( dim_num, point_num, sparse_index,
779 " SPARSE_INDEX matrix:" );
780
781 coef = new int[point_num];
782
783 webbur::sandia_sgmgg_coef_naive ( dim_num, point_num, sparse_index, coef );
784
785 webbur::i4vec_print ( point_num, coef, " COEF vector:" );
786
787 coef_sum = webbur::i4vec_sum ( point_num, coef );
788
789 std::cout << " --- -\n";
790 std::cout << " Sum: " << coef_sum << "\n";
791
792 delete [] coef;
793
794 return;
795 }
796 //****************************************************************************80
797
sandia_sgmgg_neighbor_naive_test(int dim_num,int point_num,int sparse_index[])798 void sandia_sgmgg_neighbor_naive_test ( int dim_num, int point_num,
799 int sparse_index[] )
800
801 //****************************************************************************80
802 //
803 // Purpose:
804 //
805 // SANDIA_SGMGG_NEIGHBOR_NAIVE_TEST demonstrates "naive" neighbor computations.
806 //
807 // Licensing:
808 //
809 // This code is distributed under the GNU LGPL license.
810 //
811 // Modified:
812 //
813 // 06 August 2010
814 //
815 // Author:
816 //
817 // John Burkardt
818 //
819 {
820 int *neighbor;
821
822 webbur::i4mat_transpose_print ( dim_num, point_num, sparse_index,
823 " SPARSE_INDEX matrix:" );
824
825 neighbor = new int[dim_num*point_num];
826
827 webbur::sandia_sgmgg_neighbor_naive ( dim_num, point_num, sparse_index,
828 neighbor );
829
830 webbur::i4mat_transpose_print ( dim_num, point_num, neighbor,
831 " NEIGHBOR matrix:" );
832
833 delete [] neighbor;
834
835 return;
836 }
837