1 //
2 // abstract.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <math.h>
33
34 #include <util/misc/formio.h>
35 #include <util/state/stateio.h>
36 #include <math/scmat/matrix.h>
37 #include <math/scmat/blkiter.h>
38 #include <math/scmat/elemop.h>
39
40 using namespace std;
41 using namespace sc;
42
43 /////////////////////////////////////////////////////////////////////////
44 // These member are used by the abstract SCMatrix classes.
45 /////////////////////////////////////////////////////////////////////////
46
47 /////////////////////////////////////////////////////////////////////////
48 // SCMatrixKit members
49
50 static ClassDesc SCMatrixKit_cd(
51 typeid(SCMatrixKit),"SCMatrixKit",1,"public DescribedClass",
52 0, 0, 0);
53
SCMatrixKit()54 SCMatrixKit::SCMatrixKit()
55 {
56 grp_ = MessageGrp::get_default_messagegrp();
57 }
58
SCMatrixKit(const Ref<KeyVal> & keyval)59 SCMatrixKit::SCMatrixKit(const Ref<KeyVal>& keyval)
60 {
61 grp_ << keyval->describedclassvalue("messagegrp");
62 if (grp_.null()) grp_ = MessageGrp::get_default_messagegrp();
63 }
64
~SCMatrixKit()65 SCMatrixKit::~SCMatrixKit()
66 {
67 }
68
69 SCMatrix*
restore_matrix(StateIn & s,const RefSCDimension & d1,const RefSCDimension & d2)70 SCMatrixKit::restore_matrix(StateIn& s,
71 const RefSCDimension& d1,
72 const RefSCDimension& d2)
73 {
74 SCMatrix *r = matrix(d1,d2);
75 r->restore(s);
76 return r;
77 }
78
79 SymmSCMatrix*
restore_symmmatrix(StateIn & s,const RefSCDimension & d)80 SCMatrixKit::restore_symmmatrix(StateIn& s, const RefSCDimension& d)
81 {
82 SymmSCMatrix *r = symmmatrix(d);
83 r->restore(s);
84 return r;
85 }
86
87 DiagSCMatrix*
restore_diagmatrix(StateIn & s,const RefSCDimension & d)88 SCMatrixKit::restore_diagmatrix(StateIn& s, const RefSCDimension& d)
89 {
90 DiagSCMatrix *r = diagmatrix(d);
91 r->restore(s);
92 return r;
93 }
94
95 SCVector*
restore_vector(StateIn & s,const RefSCDimension & d)96 SCMatrixKit::restore_vector(StateIn& s, const RefSCDimension& d)
97 {
98 SCVector *r = vector(d);
99 r->restore(s);
100 return r;
101 }
102
103 Ref<MessageGrp>
messagegrp() const104 SCMatrixKit::messagegrp() const
105 {
106 return grp_;
107 }
108
109 /////////////////////////////////////////////////////////////////////////
110 // SCMatrix members
111
112 static ClassDesc SCMatrix_cd(
113 typeid(SCMatrix),"SCMatrix",1,"public DescribedClass",
114 0, 0, 0);
115
SCMatrix(const RefSCDimension & a1,const RefSCDimension & a2,SCMatrixKit * kit)116 SCMatrix::SCMatrix(const RefSCDimension&a1, const RefSCDimension&a2,
117 SCMatrixKit*kit):
118 d1(a1),
119 d2(a2),
120 kit_(kit)
121 {
122 }
123
~SCMatrix()124 SCMatrix::~SCMatrix()
125 {
126 }
127
128 void
save(StateOut & s)129 SCMatrix::save(StateOut&s)
130 {
131 int nr = nrow();
132 int nc = ncol();
133 s.put(nr);
134 s.put(nc);
135 int has_subblocks = 0;
136 s.put(has_subblocks);
137 for (int i=0; i<nr; i++) {
138 for (int j=0; j<nc; j++) {
139 s.put(get_element(i,j));
140 }
141 }
142 }
143
144 void
restore(StateIn & s)145 SCMatrix::restore(StateIn& s)
146 {
147 int nrt, nr = nrow();
148 int nct, nc = ncol();
149 s.get(nrt);
150 s.get(nct);
151 if (nrt != nr || nct != nc) {
152 ExEnv::errn() << "SCMatrix::restore(): bad dimensions" << endl;
153 abort();
154 }
155 int has_subblocks;
156 s.get(has_subblocks);
157 if (!has_subblocks) {
158 for (int i=0; i<nr; i++) {
159 for (int j=0; j<nc; j++) {
160 double tmp;
161 s.get(tmp);
162 set_element(i,j, tmp);
163 }
164 }
165 }
166 else {
167 ExEnv::errn() << "SCMatrix::restore(): matrix has subblocks--cannot restore"
168 << endl;
169 abort();
170 }
171 }
172
173 double
maxabs() const174 SCMatrix::maxabs() const
175 {
176 Ref<SCElementMaxAbs> op = new SCElementMaxAbs();
177 Ref<SCElementOp> abop = op.pointer();
178 ((SCMatrix *)this)->element_op(abop);
179 return op->result();
180 }
181
182 void
randomize()183 SCMatrix::randomize()
184 {
185 Ref<SCElementOp> op = new SCElementRandomize();
186 this->element_op(op);
187 }
188
189 void
assign_val(double a)190 SCMatrix::assign_val(double a)
191 {
192 Ref<SCElementOp> op = new SCElementAssign(a);
193 this->element_op(op);
194 }
195
196 void
scale(double a)197 SCMatrix::scale(double a)
198 {
199 Ref<SCElementOp> op = new SCElementScale(a);
200 this->element_op(op);
201 }
202
203 void
scale_diagonal(double a)204 SCMatrix::scale_diagonal(double a)
205 {
206 Ref<SCElementOp> op = new SCElementScaleDiagonal(a);
207 this->element_op(op);
208 }
209
210 void
shift_diagonal(double a)211 SCMatrix::shift_diagonal(double a)
212 {
213 Ref<SCElementOp> op = new SCElementShiftDiagonal(a);
214 this->element_op(op);
215 }
216
217 void
unit()218 SCMatrix::unit()
219 {
220 this->assign(0.0);
221 this->shift_diagonal(1.0);
222 }
223
224 void
assign_r(SCMatrix * a)225 SCMatrix::assign_r(SCMatrix*a)
226 {
227 this->assign(0.0);
228 this->accumulate(a);
229 }
230
231 void
assign_p(const double * a)232 SCMatrix::assign_p(const double*a)
233 {
234 int i;
235 int nr = nrow();
236 int nc = ncol();
237 // some compilers need the following cast
238 const double **v = (const double**) new double*[nr];
239 for (i=0; i<nr; i++) {
240 v[i] = &a[i*nc];
241 }
242 assign(v);
243 delete[] v;
244 }
245
246 void
assign_pp(const double ** a)247 SCMatrix::assign_pp(const double**a)
248 {
249 int i;
250 int j;
251 int nr;
252 int nc;
253 nr = nrow();
254 nc = ncol();
255 for (i=0; i<nr; i++) {
256 for (j=0; j<nc; j++) {
257 set_element(i,j,a[i][j]);
258 }
259 }
260 }
261
262 void
convert(SCMatrix * a)263 SCMatrix::convert(SCMatrix*a)
264 {
265 assign(0.0);
266 convert_accumulate(a);
267 }
268
269 void
convert_accumulate(SCMatrix * a)270 SCMatrix::convert_accumulate(SCMatrix*a)
271 {
272 Ref<SCElementOp> op = new SCElementAccumulateSCMatrix(a);
273 element_op(op);
274 }
275
276 void
convert(double * a) const277 SCMatrix::convert(double*a) const
278 {
279 int i;
280 int nr = nrow();
281 int nc = ncol();
282 double **v = new double*[nr];
283 for (i=0; i<nr; i++) {
284 v[i] = &a[i*nc];
285 }
286 convert(v);
287 delete[] v;
288 }
289
290 void
convert(double ** a) const291 SCMatrix::convert(double**a) const
292 {
293 int i, j;
294 int nr, nc;
295 nr = nrow();
296 nc = ncol();
297 for (i=0; i<nr; i++) {
298 for (j=0; j<nc; j++) {
299 a[i][j] = get_element(i,j);
300 }
301 }
302 }
303
304 void
accumulate_product_sr(SymmSCMatrix * a,SCMatrix * b)305 SCMatrix::accumulate_product_sr(SymmSCMatrix*a,SCMatrix*b)
306 {
307 SCMatrix *t = b->copy();
308 t->transpose_this();
309 SCMatrix *t2 = this->copy();
310 t2->transpose_this();
311 t2->accumulate_product(t,a);
312 delete t;
313 t2->transpose_this();
314 assign(t2);
315 delete t2;
316 }
317
318 void
accumulate_product_dr(DiagSCMatrix * a,SCMatrix * b)319 SCMatrix::accumulate_product_dr(DiagSCMatrix*a,SCMatrix*b)
320 {
321 SCMatrix *t = b->copy();
322 t->transpose_this();
323 SCMatrix *t2 = kit()->matrix(coldim(),rowdim());
324 t2->assign(0.0);
325 t2->accumulate_product(t,a);
326 delete t;
327 t2->transpose_this();
328 accumulate(t2);
329 delete t2;
330 }
331
332 void
print(ostream & o) const333 SCMatrix::print(ostream&o) const
334 {
335 vprint(0, o, 10);
336 }
337
338 void
print(const char * t,ostream & o,int i) const339 SCMatrix::print(const char *t, ostream&o, int i) const
340 {
341 vprint(t, o, i);
342 }
343
344 SCMatrix*
clone()345 SCMatrix::clone()
346 {
347 return kit()->matrix(rowdim(),coldim());
348 }
349
350 SCMatrix*
copy()351 SCMatrix::copy()
352 {
353 SCMatrix* result = clone();
354 result->assign(this);
355 return result;
356 }
357
358 void
gen_invert_this()359 SCMatrix::gen_invert_this()
360 {
361 int i;
362
363 // Compute the singular value decomposition of this
364 RefSCMatrix U(rowdim(),rowdim(),kit());
365 RefSCMatrix V(coldim(),coldim(),kit());
366 RefSCDimension min;
367 if (coldim().n()<rowdim().n()) min = coldim();
368 else min = rowdim();
369 int nmin = min.n();
370 RefDiagSCMatrix sigma(min,kit());
371 svd_this(U.pointer(),sigma.pointer(),V.pointer());
372
373 // compute the epsilon rank of this
374 int rank = 0;
375 for (i=0; i<nmin; i++) {
376 if (fabs(sigma(i)) > 0.0000001) rank++;
377 }
378
379 RefSCDimension drank = new SCDimension(rank);
380 RefDiagSCMatrix sigma_i(drank,kit());
381 for (i=0; i<rank; i++) {
382 sigma_i(i) = 1.0/sigma(i);
383 }
384 RefSCMatrix Ur(rowdim(), drank, kit());
385 RefSCMatrix Vr(coldim(), drank, kit());
386 Ur.assign_subblock(U,0, rowdim().n()-1, 0, drank.n()-1, 0, 0);
387 Vr.assign_subblock(V,0, coldim().n()-1, 0, drank.n()-1, 0, 0);
388 assign((Vr * sigma_i * Ur.t()).t());
389 transpose_this();
390 }
391
392 void
svd_this(SCMatrix * U,DiagSCMatrix * sigma,SCMatrix * V)393 SCMatrix::svd_this(SCMatrix *U, DiagSCMatrix *sigma, SCMatrix *V)
394 {
395 ExEnv::errn() << indent << class_name() << ": SVD not implemented\n";
396 abort();
397 }
398
399 void
accumulate_product_rs(SCMatrix * a,SymmSCMatrix * b)400 SCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b)
401 {
402 SCMatrix *brect = kit()->matrix(b->dim(),b->dim());
403 brect->assign(0.0);
404 brect->accumulate(b);
405 accumulate_product(a,brect);
406 delete brect;
407 }
408
409 void
accumulate_product_ss(SymmSCMatrix * a,SymmSCMatrix * b)410 SCMatrix::accumulate_product_ss(SymmSCMatrix*a,SymmSCMatrix*b)
411 {
412 SCMatrix *arect = kit()->matrix(a->dim(),a->dim());
413 arect->assign(0.0);
414 arect->accumulate(a);
415 SCMatrix *brect = kit()->matrix(b->dim(),b->dim());
416 brect->assign(0.0);
417 brect->accumulate(b);
418 accumulate_product(arect,brect);
419 delete arect;
420 delete brect;
421 }
422
423 void
accumulate_product_rd(SCMatrix * a,DiagSCMatrix * b)424 SCMatrix::accumulate_product_rd(SCMatrix*a,DiagSCMatrix*b)
425 {
426 SCMatrix *brect = kit()->matrix(b->dim(),b->dim());
427 brect->assign(0.0);
428 brect->accumulate(b);
429 accumulate_product(a,brect);
430 delete brect;
431 }
432
433 Ref<MessageGrp>
messagegrp() const434 SCMatrix::messagegrp() const
435 {
436 return kit_->messagegrp();
437 }
438
439 /////////////////////////////////////////////////////////////////////////
440 // SymmSCMatrix member functions
441
442 static ClassDesc SymmSCMatrix_cd(
443 typeid(SymmSCMatrix),"SymmSCMatrix",1,"public DescribedClass",
444 0, 0, 0);
445
SymmSCMatrix(const RefSCDimension & a,SCMatrixKit * kit)446 SymmSCMatrix::SymmSCMatrix(const RefSCDimension&a, SCMatrixKit *kit):
447 d(a),
448 kit_(kit)
449 {
450 }
451
~SymmSCMatrix()452 SymmSCMatrix::~SymmSCMatrix()
453 {
454 }
455
456 void
save(StateOut & s)457 SymmSCMatrix::save(StateOut&s)
458 {
459 int nr = n();
460 s.put(nr);
461 for (int i=0; i<nr; i++) {
462 for (int j=0; j<=i; j++) {
463 s.put(get_element(i,j));
464 }
465 }
466 }
467
468 void
restore(StateIn & s)469 SymmSCMatrix::restore(StateIn& s)
470 {
471 int nrt, nr = n();
472 s.get(nrt);
473 if (nrt != nr) {
474 ExEnv::errn() << "SymmSCMatrix::restore(): bad dimension" << endl;
475 abort();
476 }
477 for (int i=0; i<nr; i++) {
478 for (int j=0; j<=i; j++) {
479 double tmp;
480 s.get(tmp);
481 set_element(i,j, tmp);
482 }
483 }
484 }
485
486 double
maxabs() const487 SymmSCMatrix::maxabs() const
488 {
489 Ref<SCElementMaxAbs> op = new SCElementMaxAbs();
490 Ref<SCElementOp> abop = op.pointer();
491 ((SymmSCMatrix*)this)->element_op(abop);
492 return op->result();
493 }
494
495 void
randomize()496 SymmSCMatrix::randomize()
497 {
498 Ref<SCElementOp> op = new SCElementRandomize();
499 this->element_op(op);
500 }
501
502 void
assign_val(double a)503 SymmSCMatrix::assign_val(double a)
504 {
505 Ref<SCElementOp> op = new SCElementAssign(a);
506 this->element_op(op);
507 }
508
509 void
assign_p(const double * a)510 SymmSCMatrix::assign_p(const double*a)
511 {
512 int i;
513 int nr = n();
514 // some compilers need the following cast
515 const double **v = (const double **) new double*[nr];
516 int ioff= 0;
517 for (i=0; i<nr; i++) {
518 ioff += i;
519 v[i] = &a[ioff];
520 // ioff += i;
521 }
522 assign(v);
523 delete[] v;
524 }
525
526 void
assign_pp(const double ** a)527 SymmSCMatrix::assign_pp(const double**a)
528 {
529 int i;
530 int j;
531 int nr;
532 nr = n();
533 for (i=0; i<nr; i++) {
534 for (j=0; j<=i; j++) {
535 set_element(i,j,a[i][j]);
536 }
537 }
538 }
539
540 void
convert(SymmSCMatrix * a)541 SymmSCMatrix::convert(SymmSCMatrix*a)
542 {
543 assign(0.0);
544 convert_accumulate(a);
545 }
546
547 void
convert_accumulate(SymmSCMatrix * a)548 SymmSCMatrix::convert_accumulate(SymmSCMatrix*a)
549 {
550 Ref<SCElementOp> op = new SCElementAccumulateSymmSCMatrix(a);
551 element_op(op);
552 }
553
554 void
convert(double * a) const555 SymmSCMatrix::convert(double*a) const
556 {
557 int i;
558 int nr = n();
559 double **v = new double*[nr];
560 int ioff= 0;
561 for (i=0; i<nr; i++) {
562 ioff += i;
563 v[i] = &a[ioff];
564 // ioff += i;
565 }
566 convert(v);
567 delete[] v;
568 }
569
570 void
convert(double ** a) const571 SymmSCMatrix::convert(double**a) const
572 {
573 int i;
574 int j;
575 int nr;
576 nr = n();
577 for (i=0; i<nr; i++) {
578 for (j=0; j<=i; j++) {
579 a[i][j] = get_element(i,j);
580 }
581 }
582 }
583
584 void
scale(double a)585 SymmSCMatrix::scale(double a)
586 {
587 Ref<SCElementOp> op = new SCElementScale(a);
588 this->element_op(op);
589 }
590
591 void
scale_diagonal(double a)592 SymmSCMatrix::scale_diagonal(double a)
593 {
594 Ref<SCElementOp> op = new SCElementScaleDiagonal(a);
595 this->element_op(op);
596 }
597
598 void
shift_diagonal(double a)599 SymmSCMatrix::shift_diagonal(double a)
600 {
601 Ref<SCElementOp> op = new SCElementShiftDiagonal(a);
602 this->element_op(op);
603 }
604
605 void
unit()606 SymmSCMatrix::unit()
607 {
608 this->assign(0.0);
609 this->shift_diagonal(1.0);
610 }
611
612 void
assign_s(SymmSCMatrix * a)613 SymmSCMatrix::assign_s(SymmSCMatrix*a)
614 {
615 this->assign(0.0);
616 this->accumulate(a);
617 }
618
619 void
print(ostream & o) const620 SymmSCMatrix::print(ostream&o) const
621 {
622 vprint(0, o, 10);
623 }
624
625 void
print(const char * t,ostream & o,int i) const626 SymmSCMatrix::print(const char *t, ostream&o, int i) const
627 {
628 vprint(t, o, i);
629 }
630
631 void
vprint(const char * title,ostream & out,int i) const632 SymmSCMatrix::vprint(const char* title, ostream& out, int i) const
633 {
634 RefSCMatrix m = kit()->matrix(dim(),dim());
635 m->assign(0.0);
636 m->accumulate(this);
637 m->print(title, out, i);
638 }
639
640 SymmSCMatrix*
clone()641 SymmSCMatrix::clone()
642 {
643 return kit()->symmmatrix(dim());
644 }
645
646 SymmSCMatrix*
copy()647 SymmSCMatrix::copy()
648 {
649 SymmSCMatrix* result = clone();
650 result->assign(this);
651 return result;
652 }
653
654 void
accumulate_symmetric_product(SCMatrix * a)655 SymmSCMatrix::accumulate_symmetric_product(SCMatrix *a)
656 {
657 RefSCMatrix at = a->copy();
658 at->transpose_this();
659 RefSCMatrix m = kit()->matrix(a->rowdim(),a->rowdim());
660 m->assign(0.0);
661 m->accumulate_product(a, at.pointer());
662 scale(2.0);
663 accumulate_symmetric_sum(m.pointer());
664 scale(0.5);
665 }
666
667 void
accumulate_transform(SCMatrix * a,SymmSCMatrix * b,SCMatrix::Transform t)668 SymmSCMatrix::accumulate_transform(SCMatrix *a, SymmSCMatrix *b,
669 SCMatrix::Transform t)
670 {
671 RefSCMatrix brect = kit()->matrix(b->dim(),b->dim());
672 brect->assign(0.0);
673 brect->accumulate(b);
674
675 RefSCMatrix res;
676
677 if (t == SCMatrix::TransposeTransform) {
678 RefSCMatrix at = a->copy();
679 at->transpose_this();
680
681 RefSCMatrix tmp = at->clone();
682 tmp->assign(0.0);
683
684 tmp->accumulate_product(at.pointer(), brect.pointer());
685 brect = 0;
686 at = 0;
687
688 res = kit()->matrix(dim(),dim());
689 res->assign(0.0);
690 res->accumulate_product(tmp.pointer(), a);
691 }
692 else {
693 RefSCMatrix tmp = a->clone();
694 tmp->assign(0.0);
695
696 tmp->accumulate_product(a,brect);
697 brect = 0;
698
699 RefSCMatrix at = a->copy();
700 at->transpose_this();
701
702 res = kit()->matrix(dim(),dim());
703 res->assign(0.0);
704 res->accumulate_product(tmp.pointer(), at.pointer());
705 at = 0;
706 }
707
708 scale(2.0);
709 accumulate_symmetric_sum(res.pointer());
710 scale(0.5);
711 }
712
713 void
accumulate_transform(SCMatrix * a,DiagSCMatrix * b,SCMatrix::Transform t)714 SymmSCMatrix::accumulate_transform(SCMatrix *a, DiagSCMatrix *b,
715 SCMatrix::Transform t)
716 {
717 RefSCMatrix m = kit()->matrix(a->rowdim(),a->rowdim());
718 RefSCMatrix brect = kit()->matrix(b->dim(),b->dim());
719 brect->assign(0.0);
720 brect->accumulate(b);
721
722 RefSCMatrix tmp = a->clone();
723 tmp->assign(0.0);
724
725 RefSCMatrix res;
726
727 if (t == SCMatrix::TransposeTransform) {
728 RefSCMatrix at = a->copy();
729 at->transpose_this();
730
731 tmp->accumulate_product(at.pointer(), brect.pointer());
732 brect = 0;
733 at = 0;
734
735 res = kit()->matrix(dim(),dim());
736 res->assign(0.0);
737 res->accumulate_product(tmp.pointer(), a);
738 }
739 else {
740 tmp->accumulate_product(a, brect.pointer());
741 brect = 0;
742
743 RefSCMatrix at = a->copy();
744 at->transpose_this();
745
746 res = kit()->matrix(dim(),dim());
747 res->assign(0.0);
748 res->accumulate_product(tmp.pointer(), at.pointer());
749 at = 0;
750 }
751
752 tmp = 0;
753
754 scale(2.0);
755 accumulate_symmetric_sum(res.pointer());
756 scale(0.5);
757 }
758
759 void
accumulate_transform(SymmSCMatrix * a,SymmSCMatrix * b)760 SymmSCMatrix::accumulate_transform(SymmSCMatrix *a, SymmSCMatrix *b)
761 {
762 RefSCMatrix m = kit()->matrix(a->dim(),a->dim());
763 m->assign(0.0);
764 m->accumulate(a);
765 accumulate_transform(m.pointer(),b);
766 }
767
768 void
accumulate_symmetric_outer_product(SCVector * v)769 SymmSCMatrix::accumulate_symmetric_outer_product(SCVector *v)
770 {
771 RefSCMatrix m = kit()->matrix(dim(),dim());
772 m->assign(0.0);
773 m->accumulate_outer_product(v,v);
774
775 scale(2.0);
776 accumulate_symmetric_sum(m.pointer());
777 scale(0.5);
778 }
779
780 double
scalar_product(SCVector * v)781 SymmSCMatrix::scalar_product(SCVector *v)
782 {
783 RefSCVector v2 = kit()->vector(dim());
784 v2->assign(0.0);
785 v2->accumulate_product(this,v);
786 return v2->scalar_product(v);
787 }
788
789 Ref<MessageGrp>
messagegrp() const790 SymmSCMatrix::messagegrp() const
791 {
792 return kit_->messagegrp();
793 }
794
795 /////////////////////////////////////////////////////////////////////////
796 // DiagSCMatrix member functions
797
798 static ClassDesc DiagSCMatrix_cd(
799 typeid(DiagSCMatrix),"DiagSCMatrix",1,"public DescribedClass",
800 0, 0, 0);
801
DiagSCMatrix(const RefSCDimension & a,SCMatrixKit * kit)802 DiagSCMatrix::DiagSCMatrix(const RefSCDimension&a, SCMatrixKit *kit):
803 d(a),
804 kit_(kit)
805 {
806 }
807
~DiagSCMatrix()808 DiagSCMatrix::~DiagSCMatrix()
809 {
810 }
811
812 void
save(StateOut & s)813 DiagSCMatrix::save(StateOut&s)
814 {
815 int nr = n();
816 s.put(nr);
817 for (int i=0; i<nr; i++) {
818 s.put(get_element(i));
819 }
820 }
821
822 void
restore(StateIn & s)823 DiagSCMatrix::restore(StateIn& s)
824 {
825 int nrt, nr = n();
826 s.get(nrt);
827 if (nrt != nr) {
828 ExEnv::errn() << "DiagSCMatrix::restore(): bad dimension" << endl;
829 abort();
830 }
831 for (int i=0; i<nr; i++) {
832 double tmp;
833 s.get(tmp);
834 set_element(i, tmp);
835 }
836 }
837
838 double
maxabs() const839 DiagSCMatrix::maxabs() const
840 {
841 Ref<SCElementMaxAbs> op = new SCElementMaxAbs();
842 Ref<SCElementOp> abop = op.pointer();
843 ((DiagSCMatrix*)this)->element_op(abop);
844 return op->result();
845 }
846
847 void
randomize()848 DiagSCMatrix::randomize()
849 {
850 Ref<SCElementOp> op = new SCElementRandomize();
851 this->element_op(op);
852 }
853
854 void
assign_val(double a)855 DiagSCMatrix::assign_val(double a)
856 {
857 Ref<SCElementOp> op = new SCElementAssign(a);
858 this->element_op(op);
859 }
860
861 void
assign_p(const double * a)862 DiagSCMatrix::assign_p(const double*a)
863 {
864 int i;
865 int nr = n();
866 for (i=0; i<nr; i++) {
867 set_element(i,a[i]);
868 }
869 }
870
871 void
convert(DiagSCMatrix * a)872 DiagSCMatrix::convert(DiagSCMatrix*a)
873 {
874 assign(0.0);
875 convert_accumulate(a);
876 }
877
878 void
convert_accumulate(DiagSCMatrix * a)879 DiagSCMatrix::convert_accumulate(DiagSCMatrix*a)
880 {
881 Ref<SCElementOp> op = new SCElementAccumulateDiagSCMatrix(a);
882 element_op(op);
883 }
884
885 void
convert(double * a) const886 DiagSCMatrix::convert(double*a) const
887 {
888 int i;
889 int nr = n();
890 for (i=0; i<nr; i++) {
891 a[i] = get_element(i);
892 }
893 }
894
895 void
scale(double a)896 DiagSCMatrix::scale(double a)
897 {
898 Ref<SCElementOp> op = new SCElementScale(a);
899 this->element_op(op);
900 }
901
902 void
assign_d(DiagSCMatrix * a)903 DiagSCMatrix::assign_d(DiagSCMatrix*a)
904 {
905 this->assign(0.0);
906 this->accumulate(a);
907 }
908
909 void
print(ostream & o) const910 DiagSCMatrix::print(ostream&o) const
911 {
912 vprint(0, o, 10);
913 }
914
915 void
print(const char * t,ostream & o,int i) const916 DiagSCMatrix::print(const char *t, ostream&o, int i) const
917 {
918 vprint(t, o, i);
919 }
920
921 void
vprint(const char * title,ostream & out,int i) const922 DiagSCMatrix::vprint(const char* title, ostream& out, int i) const
923 {
924 RefSCMatrix m = kit()->matrix(dim(),dim());
925 m->assign(0.0);
926 m->accumulate(this);
927 m->print(title, out, i);
928 }
929
930 DiagSCMatrix*
clone()931 DiagSCMatrix::clone()
932 {
933 return kit()->diagmatrix(dim());
934 }
935
936 DiagSCMatrix*
copy()937 DiagSCMatrix::copy()
938 {
939 DiagSCMatrix* result = clone();
940 result->assign(this);
941 return result;
942 }
943
944 Ref<MessageGrp>
messagegrp() const945 DiagSCMatrix::messagegrp() const
946 {
947 return kit_->messagegrp();
948 }
949
950 /////////////////////////////////////////////////////////////////////////
951 // These member are used by the abstract SCVector classes.
952 /////////////////////////////////////////////////////////////////////////
953
954 static ClassDesc SCVector_cd(
955 typeid(SCVector),"SCVector",1,"public DescribedClass",
956 0, 0, 0);
957
SCVector(const RefSCDimension & a,SCMatrixKit * kit)958 SCVector::SCVector(const RefSCDimension&a, SCMatrixKit *kit):
959 d(a),
960 kit_(kit)
961 {
962 }
963
~SCVector()964 SCVector::~SCVector()
965 {
966 }
967
968 void
save(StateOut & s)969 SCVector::save(StateOut&s)
970 {
971 int nr = n();
972 s.put(nr);
973 for (int i=0; i<nr; i++) {
974 s.put(get_element(i));
975 }
976 }
977
978 void
restore(StateIn & s)979 SCVector::restore(StateIn& s)
980 {
981 int nrt, nr = n();
982 s.get(nrt);
983 if (nrt != nr) {
984 ExEnv::errn() << "SCVector::restore(): bad dimension" << endl;
985 abort();
986 }
987 for (int i=0; i<nr; i++) {
988 double tmp;
989 s.get(tmp);
990 set_element(i, tmp);
991 }
992 }
993
994 double
maxabs() const995 SCVector::maxabs() const
996 {
997 Ref<SCElementMaxAbs> op = new SCElementMaxAbs();
998 Ref<SCElementOp> abop = op.pointer();
999 ((SCVector*)this)->element_op(abop);
1000 return op->result();
1001 }
1002
1003 void
randomize()1004 SCVector::randomize()
1005 {
1006 Ref<SCElementOp> op = new SCElementRandomize();
1007 this->element_op(op);
1008 }
1009
1010 void
assign_val(double a)1011 SCVector::assign_val(double a)
1012 {
1013 Ref<SCElementOp> op = new SCElementAssign(a);
1014 this->element_op(op);
1015 }
1016
1017 void
assign_p(const double * a)1018 SCVector::assign_p(const double*a)
1019 {
1020 int i;
1021 int nr = n();
1022 for (i=0; i<nr; i++) {
1023 set_element(i,a[i]);
1024 }
1025 }
1026
1027 void
convert(SCVector * a)1028 SCVector::convert(SCVector*a)
1029 {
1030 assign(0.0);
1031 convert_accumulate(a);
1032 }
1033
1034 void
convert_accumulate(SCVector * a)1035 SCVector::convert_accumulate(SCVector*a)
1036 {
1037 Ref<SCElementOp> op = new SCElementAccumulateSCVector(a);
1038 element_op(op);
1039 }
1040
1041 void
convert(double * a) const1042 SCVector::convert(double*a) const
1043 {
1044 int i;
1045 int nr = n();
1046 for (i=0; i<nr; i++) {
1047 a[i] = get_element(i);
1048 }
1049 }
1050
1051 void
scale(double a)1052 SCVector::scale(double a)
1053 {
1054 Ref<SCElementOp> op = new SCElementScale(a);
1055 this->element_op(op);
1056 }
1057
1058 void
assign_v(SCVector * a)1059 SCVector::assign_v(SCVector*a)
1060 {
1061 this->assign(0.0);
1062 this->accumulate(a);
1063 }
1064
1065 void
print(ostream & o) const1066 SCVector::print(ostream&o) const
1067 {
1068 vprint(0, o, 10);
1069 }
1070
1071 void
print(const char * t,ostream & o,int i) const1072 SCVector::print(const char *t, ostream&o, int i) const
1073 {
1074 vprint(t, o, i);
1075 }
1076
1077 void
normalize()1078 SCVector::normalize()
1079 {
1080 double norm = sqrt(scalar_product(this));
1081 if (norm > 1.e-20) norm = 1.0/norm;
1082 else {
1083 ExEnv::errn() << indent
1084 << "SCVector::normalize: tried to normalize tiny vector\n";
1085 abort();
1086 }
1087 scale(norm);
1088 }
1089
1090 SCVector*
clone()1091 SCVector::clone()
1092 {
1093 return kit()->vector(dim());
1094 }
1095
1096 SCVector*
copy()1097 SCVector::copy()
1098 {
1099 SCVector* result = clone();
1100 result->assign(this);
1101 return result;
1102 }
1103
1104 void
accumulate_product_sv(SymmSCMatrix * m,SCVector * v)1105 SCVector::accumulate_product_sv(SymmSCMatrix *m, SCVector *v)
1106 {
1107 RefSCMatrix mrect = kit()->matrix(m->dim(),m->dim());
1108 mrect->assign(0.0);
1109 mrect->accumulate(m);
1110 accumulate_product(mrect.pointer(), v);
1111 }
1112
1113 Ref<MessageGrp>
messagegrp() const1114 SCVector::messagegrp() const
1115 {
1116 return kit_->messagegrp();
1117 }
1118
1119 /////////////////////////////////////////////////////////////////////////////
1120
1121 // Local Variables:
1122 // mode: c++
1123 // c-file-style: "CLJ"
1124 // End:
1125