Initial revision
[blender.git] / intern / iksolver / intern / TNT / triang.h
1 /**
2  * $Id$
3  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version. The Blender
9  * Foundation also sells licenses for use in proprietary software under
10  * the Blender License.  See http://www.blender.org/BL/ for information
11  * about this.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
21  *
22  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): none yet.
28  *
29  * ***** END GPL/BL DUAL LICENSE BLOCK *****
30  */
31
32 /*
33
34 *
35 * Template Numerical Toolkit (TNT): Linear Algebra Module
36 *
37 * Mathematical and Computational Sciences Division
38 * National Institute of Technology,
39 * Gaithersburg, MD USA
40 *
41 *
42 * This software was developed at the National Institute of Standards and
43 * Technology (NIST) by employees of the Federal Government in the course
44 * of their official duties. Pursuant to title 17 Section 105 of the
45 * United States Code, this software is not subject to copyright protection
46 * and is in the public domain.  The Template Numerical Toolkit (TNT) is
47 * an experimental system.  NIST assumes no responsibility whatsoever for
48 * its use by other parties, and makes no guarantees, expressed or implied,
49 * about its quality, reliability, or any other characteristic.
50 *
51 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
52 * see http://math.nist.gov/tnt for latest updates.
53 *
54 */
55
56
57
58 // Triangular Matrices (Views and Adpators)
59
60 #ifndef TRIANG_H
61 #define TRIANG_H
62
63 // default to use lower-triangular portions of arrays
64 // for symmetric matrices.
65
66 namespace TNT
67 {
68
69 template <class MaTRiX>
70 class LowerTriangularView
71 {
72     protected:
73
74
75         const MaTRiX  &A_;
76         const typename MaTRiX::element_type zero_;
77
78     public:
79
80
81     typedef typename MaTRiX::const_reference const_reference;
82     typedef const typename MaTRiX::element_type element_type;
83     typedef const typename MaTRiX::element_type value_type;
84     typedef element_type T;
85
86     Subscript dim(Subscript d) const {  return A_.dim(d); }
87     Subscript lbound() const { return A_.lbound(); }
88     Subscript num_rows() const { return A_.num_rows(); }
89     Subscript num_cols() const { return A_.num_cols(); }
90     
91     
92     // constructors
93
94     LowerTriangularView(/*const*/ MaTRiX &A) : A_(A),  zero_(0) {}
95
96
97     inline const_reference get(Subscript i, Subscript j) const
98     { 
99 #ifdef TNT_BOUNDS_CHECK
100         assert(lbound()<=i);
101         assert(i<=A_.num_rows() + lbound() - 1);
102         assert(lbound()<=j);
103         assert(j<=A_.num_cols() + lbound() - 1);
104 #endif
105         if (i<j) 
106             return zero_;
107         else
108             return A_(i,j);
109     }
110
111
112     inline const_reference operator() (Subscript i, Subscript j) const
113     {
114 #ifdef TNT_BOUNDS_CHECK
115         assert(lbound()<=i);
116         assert(i<=A_.num_rows() + lbound() - 1);
117         assert(lbound()<=j);
118         assert(j<=A_.num_cols() + lbound() - 1);
119 #endif
120         if (i<j) 
121             return zero_;
122         else
123             return A_(i,j);
124     }
125
126 #ifdef TNT_USE_REGIONS 
127
128     typedef const_Region2D< LowerTriangularView<MaTRiX> > 
129                     const_Region;
130
131     const_Region operator()(/*const*/ Index1D &I,
132             /*const*/ Index1D &J) const
133     {
134         return const_Region(*this, I, J);
135     }
136
137     const_Region operator()(Subscript i1, Subscript i2,
138             Subscript j1, Subscript j2) const
139     {
140         return const_Region(*this, i1, i2, j1, j2);
141     }
142
143
144
145 #endif
146 // TNT_USE_REGIONS
147
148 };
149
150
151 /* *********** Lower_triangular_view() algorithms ****************** */
152
153 template <class MaTRiX, class VecToR>
154 VecToR matmult(/*const*/ LowerTriangularView<MaTRiX> &A, VecToR &x)
155 {
156     Subscript M = A.num_rows();
157     Subscript N = A.num_cols();
158
159     assert(N == x.dim());
160
161     Subscript i, j;
162     typename MaTRiX::element_type sum=0.0;
163     VecToR result(M);
164
165     Subscript start = A.lbound();
166     Subscript Mend = M + A.lbound() -1 ;
167
168     for (i=start; i<=Mend; i++)
169     {
170         sum = 0.0;
171         for (j=start; j<=i; j++)
172             sum = sum + A(i,j)*x(j);
173         result(i) = sum;
174     }
175
176     return result;
177 }
178
179 template <class MaTRiX, class VecToR>
180 inline VecToR operator*(/*const*/ LowerTriangularView<MaTRiX> &A, VecToR &x)
181 {
182     return matmult(A,x);
183 }
184
185 template <class MaTRiX>
186 class UnitLowerTriangularView
187 {
188     protected:
189
190         const MaTRiX  &A_;
191         const typename MaTRiX::element_type zero;
192         const typename MaTRiX::element_type one;
193
194     public:
195
196     typedef typename MaTRiX::const_reference const_reference;
197     typedef typename MaTRiX::element_type element_type;
198     typedef typename MaTRiX::element_type value_type;
199     typedef element_type T;
200
201     Subscript lbound() const { return 1; }
202     Subscript dim(Subscript d) const {  return A_.dim(d); }
203     Subscript num_rows() const { return A_.num_rows(); }
204     Subscript num_cols() const { return A_.num_cols(); }
205
206     
207     // constructors
208
209     UnitLowerTriangularView(/*const*/ MaTRiX &A) : A_(A), zero(0), one(1) {}
210
211
212     inline const_reference get(Subscript i, Subscript j) const
213     { 
214 #ifdef TNT_BOUNDS_CHECK
215         assert(1<=i);
216         assert(i<=A_.dim(1));
217         assert(1<=j);
218         assert(j<=A_.dim(2));
219         assert(0<=i && i<A_.dim(0) && 0<=j && j<A_.dim(1));
220 #endif
221         if (i>j)
222             return A_(i,j);
223         else if (i==j)
224             return one;
225         else 
226             return zero;
227     }
228
229
230     inline const_reference operator() (Subscript i, Subscript j) const
231     {
232 #ifdef TNT_BOUNDS_CHECK
233         assert(1<=i);
234         assert(i<=A_.dim(1));
235         assert(1<=j);
236         assert(j<=A_.dim(2));
237 #endif
238         if (i>j)
239             return A_(i,j);
240         else if (i==j)
241             return one;
242         else 
243             return zero;
244     }
245
246
247 #ifdef TNT_USE_REGIONS 
248   // These are the "index-aware" features
249
250     typedef const_Region2D< UnitLowerTriangularView<MaTRiX> > 
251                     const_Region;
252
253     const_Region operator()(/*const*/ Index1D &I,
254             /*const*/ Index1D &J) const
255     {
256         return const_Region(*this, I, J);
257     }
258
259     const_Region operator()(Subscript i1, Subscript i2,
260             Subscript j1, Subscript j2) const
261     {
262         return const_Region(*this, i1, i2, j1, j2);
263     }
264 #endif
265 // TNT_USE_REGIONS
266 };
267
268 template <class MaTRiX>
269 LowerTriangularView<MaTRiX> Lower_triangular_view(
270     /*const*/ MaTRiX &A)
271 {
272     return LowerTriangularView<MaTRiX>(A);
273 }
274
275
276 template <class MaTRiX>
277 UnitLowerTriangularView<MaTRiX> Unit_lower_triangular_view(
278     /*const*/ MaTRiX &A)
279 {
280     return UnitLowerTriangularView<MaTRiX>(A);
281 }
282
283 template <class MaTRiX, class VecToR>
284 VecToR matmult(/*const*/ UnitLowerTriangularView<MaTRiX> &A, VecToR &x)
285 {
286     Subscript M = A.num_rows();
287     Subscript N = A.num_cols();
288
289     assert(N == x.dim());
290
291     Subscript i, j;
292     typename MaTRiX::element_type sum=0.0;
293     VecToR result(M);
294
295     Subscript start = A.lbound();
296     Subscript Mend = M + A.lbound() -1 ;
297
298     for (i=start; i<=Mend; i++)
299     {
300         sum = 0.0;
301         for (j=start; j<i; j++)
302             sum = sum + A(i,j)*x(j);
303         result(i) = sum + x(i);
304     }
305
306     return result;
307 }
308
309 template <class MaTRiX, class VecToR>
310 inline VecToR operator*(/*const*/ UnitLowerTriangularView<MaTRiX> &A, VecToR &x)
311 {
312     return matmult(A,x);
313 }
314
315
316 //********************** Algorithms *************************************
317
318
319
320 template <class MaTRiX>
321 std::ostream& operator<<(std::ostream &s, const LowerTriangularView<MaTRiX>&A)
322 {
323     Subscript M=A.num_rows();
324     Subscript N=A.num_cols();
325
326     s << M << " " << N << endl;
327
328     for (Subscript i=1; i<=M; i++)
329     {
330         for (Subscript j=1; j<=N; j++)
331         {
332             s << A(i,j) << " ";
333         }
334         s << endl;
335     }
336
337
338     return s;
339 }
340
341 template <class MaTRiX>
342 std::ostream& operator<<(std::ostream &s, 
343     const UnitLowerTriangularView<MaTRiX>&A)
344 {
345     Subscript M=A.num_rows();
346     Subscript N=A.num_cols();
347
348     s << M << " " << N << endl;
349
350     for (Subscript i=1; i<=M; i++)
351     {
352         for (Subscript j=1; j<=N; j++)
353         {
354             s << A(i,j) << " ";
355         }
356         s << endl;
357     }
358
359
360     return s;
361 }
362
363
364
365 // ******************* Upper Triangular Section **************************
366
367 template <class MaTRiX>
368 class UpperTriangularView
369 {
370     protected:
371
372
373         /*const*/ MaTRiX  &A_;
374         /*const*/ typename MaTRiX::element_type zero_;
375
376     public:
377
378
379     typedef typename MaTRiX::const_reference const_reference;
380     typedef /*const*/ typename MaTRiX::element_type element_type;
381     typedef /*const*/ typename MaTRiX::element_type value_type;
382     typedef element_type T;
383
384     Subscript dim(Subscript d) const {  return A_.dim(d); }
385     Subscript lbound() const { return A_.lbound(); }
386     Subscript num_rows() const { return A_.num_rows(); }
387     Subscript num_cols() const { return A_.num_cols(); }
388     
389     
390     // constructors
391
392     UpperTriangularView(/*const*/ MaTRiX &A) : A_(A),  zero_(0) {}
393
394
395     inline const_reference get(Subscript i, Subscript j) const
396     { 
397 #ifdef TNT_BOUNDS_CHECK
398         assert(lbound()<=i);
399         assert(i<=A_.num_rows() + lbound() - 1);
400         assert(lbound()<=j);
401         assert(j<=A_.num_cols() + lbound() - 1);
402 #endif
403         if (i>j) 
404             return zero_;
405         else
406             return A_(i,j);
407     }
408
409
410     inline const_reference operator() (Subscript i, Subscript j) const
411     {
412 #ifdef TNT_BOUNDS_CHECK
413         assert(lbound()<=i);
414         assert(i<=A_.num_rows() + lbound() - 1);
415         assert(lbound()<=j);
416         assert(j<=A_.num_cols() + lbound() - 1);
417 #endif
418         if (i>j) 
419             return zero_;
420         else
421             return A_(i,j);
422     }
423
424 #ifdef TNT_USE_REGIONS 
425
426     typedef const_Region2D< UpperTriangularView<MaTRiX> > 
427                     const_Region;
428
429     const_Region operator()(const Index1D &I,
430             const Index1D &J) const
431     {
432         return const_Region(*this, I, J);
433     }
434
435     const_Region operator()(Subscript i1, Subscript i2,
436             Subscript j1, Subscript j2) const
437     {
438         return const_Region(*this, i1, i2, j1, j2);
439     }
440
441
442
443 #endif
444 // TNT_USE_REGIONS
445
446 };
447
448
449 /* *********** Upper_triangular_view() algorithms ****************** */
450
451 template <class MaTRiX, class VecToR>
452 VecToR matmult(/*const*/ UpperTriangularView<MaTRiX> &A, VecToR &x)
453 {
454     Subscript M = A.num_rows();
455     Subscript N = A.num_cols();
456
457     assert(N == x.dim());
458
459     Subscript i, j;
460     typename VecToR::element_type sum=0.0;
461     VecToR result(M);
462
463     Subscript start = A.lbound();
464     Subscript Mend = M + A.lbound() -1 ;
465
466     for (i=start; i<=Mend; i++)
467     {
468         sum = 0.0;
469         for (j=i; j<=N; j++)
470             sum = sum + A(i,j)*x(j);
471         result(i) = sum;
472     }
473
474     return result;
475 }
476
477 template <class MaTRiX, class VecToR>
478 inline VecToR operator*(/*const*/ UpperTriangularView<MaTRiX> &A, VecToR &x)
479 {
480     return matmult(A,x);
481 }
482
483 template <class MaTRiX>
484 class UnitUpperTriangularView
485 {
486     protected:
487
488         const MaTRiX  &A_;
489         const typename MaTRiX::element_type zero;
490         const typename MaTRiX::element_type one;
491
492     public:
493
494     typedef typename MaTRiX::const_reference const_reference;
495     typedef typename MaTRiX::element_type element_type;
496     typedef typename MaTRiX::element_type value_type;
497     typedef element_type T;
498
499     Subscript lbound() const { return 1; }
500     Subscript dim(Subscript d) const {  return A_.dim(d); }
501     Subscript num_rows() const { return A_.num_rows(); }
502     Subscript num_cols() const { return A_.num_cols(); }
503
504     
505     // constructors
506
507     UnitUpperTriangularView(/*const*/ MaTRiX &A) : A_(A), zero(0), one(1) {}
508
509
510     inline const_reference get(Subscript i, Subscript j) const
511     { 
512 #ifdef TNT_BOUNDS_CHECK
513         assert(1<=i);
514         assert(i<=A_.dim(1));
515         assert(1<=j);
516         assert(j<=A_.dim(2));
517         assert(0<=i && i<A_.dim(0) && 0<=j && j<A_.dim(1));
518 #endif
519         if (i<j)
520             return A_(i,j);
521         else if (i==j)
522             return one;
523         else 
524             return zero;
525     }
526
527
528     inline const_reference operator() (Subscript i, Subscript j) const
529     {
530 #ifdef TNT_BOUNDS_CHECK
531         assert(1<=i);
532         assert(i<=A_.dim(1));
533         assert(1<=j);
534         assert(j<=A_.dim(2));
535 #endif
536         if (i<j)
537             return A_(i,j);
538         else if (i==j)
539             return one;
540         else 
541             return zero;
542     }
543
544
545 #ifdef TNT_USE_REGIONS 
546   // These are the "index-aware" features
547
548     typedef const_Region2D< UnitUpperTriangularView<MaTRiX> > 
549                     const_Region;
550
551     const_Region operator()(const Index1D &I,
552             const Index1D &J) const
553     {
554         return const_Region(*this, I, J);
555     }
556
557     const_Region operator()(Subscript i1, Subscript i2,
558             Subscript j1, Subscript j2) const
559     {
560         return const_Region(*this, i1, i2, j1, j2);
561     }
562 #endif
563 // TNT_USE_REGIONS
564 };
565
566 template <class MaTRiX>
567 UpperTriangularView<MaTRiX> Upper_triangular_view(
568     /*const*/ MaTRiX &A)
569 {
570     return UpperTriangularView<MaTRiX>(A);
571 }
572
573
574 template <class MaTRiX>
575 UnitUpperTriangularView<MaTRiX> Unit_upper_triangular_view(
576     /*const*/ MaTRiX &A)
577 {
578     return UnitUpperTriangularView<MaTRiX>(A);
579 }
580
581 template <class MaTRiX, class VecToR>
582 VecToR matmult(/*const*/ UnitUpperTriangularView<MaTRiX> &A, VecToR &x)
583 {
584     Subscript M = A.num_rows();
585     Subscript N = A.num_cols();
586
587     assert(N == x.dim());
588
589     Subscript i, j;
590     typename VecToR::element_type sum=0.0;
591     VecToR result(M);
592
593     Subscript start = A.lbound();
594     Subscript Mend = M + A.lbound() -1 ;
595
596     for (i=start; i<=Mend; i++)
597     {
598         sum = x(i);
599         for (j=i+1; j<=N; j++)
600             sum = sum + A(i,j)*x(j);
601         result(i) = sum + x(i);
602     }
603
604     return result;
605 }
606
607 template <class MaTRiX, class VecToR>
608 inline VecToR operator*(/*const*/ UnitUpperTriangularView<MaTRiX> &A, VecToR &x)
609 {
610     return matmult(A,x);
611 }
612
613
614 //********************** Algorithms *************************************
615
616
617
618 template <class MaTRiX>
619 std::ostream& operator<<(std::ostream &s, 
620     /*const*/ UpperTriangularView<MaTRiX>&A)
621 {
622     Subscript M=A.num_rows();
623     Subscript N=A.num_cols();
624
625     s << M << " " << N << endl;
626
627     for (Subscript i=1; i<=M; i++)
628     {
629         for (Subscript j=1; j<=N; j++)
630         {
631             s << A(i,j) << " ";
632         }
633         s << endl;
634     }
635
636
637     return s;
638 }
639
640 template <class MaTRiX>
641 std::ostream& operator<<(std::ostream &s, 
642         /*const*/ UnitUpperTriangularView<MaTRiX>&A)
643 {
644     Subscript M=A.num_rows();
645     Subscript N=A.num_cols();
646
647     s << M << " " << N << endl;
648
649     for (Subscript i=1; i<=M; i++)
650     {
651         for (Subscript j=1; j<=N; j++)
652         {
653             s << A(i,j) << " ";
654         }
655         s << endl;
656     }
657
658
659     return s;
660 }
661
662 } // namespace TNT
663
664
665
666
667
668 #endif 
669 //TRIANG_H
670