source: ThirdParty/mpqc_open/src/lib/math/scmat/matrix.cc@ 13e5be

stable v1.7.0
Last change on this file since 13e5be was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 32.7 KB
RevLine 
[0b990d]1//
2// matrix.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 <util/state/stateio.h>
33#include <math/scmat/block.h>
34#include <math/scmat/blkiter.h>
35#include <math/scmat/matrix.h>
36#include <math/scmat/blocked.h>
37
38using namespace std;
39
40namespace sc {
41
42/////////////////////////////////////////////////////////////////////////////
43// SCDimension reference member functions
44
45RefSCDimension::RefSCDimension() {}
46
47RefSCDimension::RefSCDimension (const RefSCDimension & o):
48 Ref<SCDimension> (o) {}
49
50RefSCDimension::RefSCDimension (SCDimension * o): Ref<SCDimension> (o) {}
51
52RefSCDimension::~RefSCDimension () {}
53
54RefSCDimension&
55RefSCDimension::operator=(SCDimension* cr)
56{
57 Ref<SCDimension>::operator=(cr);
58 return *this;
59}
60
61RefSCDimension&
62RefSCDimension::operator<<(const RefBase & c)
63{
64 Ref<SCDimension>::operator<<(c);
65 return *this;
66}
67
68RefSCDimension&
69RefSCDimension::operator<<(RefCount*a)
70{
71 Ref<SCDimension>::operator<<(a);
72 return *this;
73}
74
75RefSCDimension&
76RefSCDimension::operator=(const RefSCDimension & c)
77{
78 Ref<SCDimension>::operator=(c);
79 return *this;
80}
81
82int
83RefSCDimension::n() const
84{
85 int result;
86 if (null()) result = 0;
87 else result = pointer()->n();
88 return result;
89}
90
91RefSCDimension::operator int() const
92{
93 if (null()) return 0;
94 return pointer()->n();
95}
96
97/////////////////////////////////////////////////////////////////////////////
98// SCMatrix reference member functions
99RefSCMatrix::RefSCMatrix() {}
100
101RefSCMatrix::RefSCMatrix (const RefSCMatrix & o): Ref<SCMatrix> (o) {}
102
103RefSCMatrix::RefSCMatrix (SCMatrix * o): Ref<SCMatrix> (o) {}
104
105RefSCMatrix::~RefSCMatrix () {}
106
107RefSCMatrix&
108RefSCMatrix::operator=(SCMatrix* cr)
109{
110 Ref<SCMatrix>::operator=(cr);
111 return *this;
112}
113
114RefSCMatrix&
115RefSCMatrix::operator=(const RefSCMatrix & c)
116{
117 Ref<SCMatrix>::operator=(c);
118 return *this;
119}
120
121RefSCMatrix::RefSCMatrix(const RefSCDimension&a,const RefSCDimension&b,
122 const Ref<SCMatrixKit>&k)
123{
124 assign_pointer(k->matrix(a,b));
125}
126
127void
128RefSCMatrix::set_element(int i, int j, double a) const
129{
130 require_nonnull();
131 pointer()->set_element(i,j,a);
132}
133
134void
135RefSCMatrix::accumulate_element(int i, int j, double a) const
136{
137 require_nonnull();
138 pointer()->accumulate_element(i,j,a);
139}
140
141double
142RefSCMatrix::get_element(int i, int j) const
143{
144 require_nonnull();
145 return pointer()->get_element(i,j);
146}
147
148RefSCVector
149RefSCMatrix::operator*(const RefSCVector&a) const
150{
151 require_nonnull();
152 a.require_nonnull();
153
154 RefSCVector r = kit()->vector(rowdim());
155 r->assign(0.0);
156 r->accumulate_product(pointer(),a.pointer());
157 return r;
158}
159
160RefSCMatrix
161RefSCMatrix::operator*(const RefSCMatrix&a) const
162{
163 require_nonnull();
164 a.require_nonnull();
165
166 RefSCMatrix r = kit()->matrix(rowdim(),a->coldim());
167 r->assign(0.0);
168 r->accumulate_product(pointer(),a.pointer());
169 return r;
170}
171
172RefSCMatrix
173RefSCMatrix::operator*(const RefSymmSCMatrix&a) const
174{
175 require_nonnull();
176 a.require_nonnull();
177
178 RefSCMatrix r = kit()->matrix(rowdim(),a->dim());
179 r->assign(0.0);
180 r->accumulate_product(pointer(),a.pointer());
181 return r;
182}
183
184RefSCMatrix
185RefSCMatrix::operator*(const RefDiagSCMatrix&a) const
186{
187 require_nonnull();
188 a.require_nonnull();
189
190 RefSCMatrix r = kit()->matrix(rowdim(),a->dim());
191 r->assign(0.0);
192 r->accumulate_product(pointer(),a.pointer());
193 return r;
194}
195
196RefSCMatrix
197RefSCMatrix::operator+(const RefSCMatrix&a) const
198{
199 require_nonnull();
200 a.require_nonnull();
201
202 RefSCMatrix ret(rowdim(),coldim(),kit());
203
204 ret->assign(pointer());
205 ret->accumulate(a.pointer());
206
207 return ret;
208}
209
210RefSCMatrix
211RefSCMatrix::operator-(const RefSCMatrix&a) const
212{
213 require_nonnull();
214 a.require_nonnull();
215
216 RefSCMatrix ret(rowdim(),coldim(),kit());
217
218 ret->assign(a.pointer());
219 ret->scale(-1.0);
220 ret->accumulate(pointer());
221
222 return ret;
223}
224
225RefSCMatrix
226RefSCMatrix::t() const
227{
228 require_nonnull();
229
230 RefSCMatrix ret;
231 ret = clone();
232 ret->assign(pointer());
233 ret->transpose_this();
234 return ret;
235}
236
237RefSCMatrix
238RefSCMatrix::i() const
239{
240 require_nonnull();
241
242 RefSCMatrix ret;
243 ret = clone();
244 ret->assign(pointer());
245 ret->invert_this();
246 return ret;
247}
248
249RefSCMatrix
250RefSCMatrix::gi() const
251{
252 require_nonnull();
253
254 RefSCMatrix ret;
255 ret = clone();
256 ret->assign(pointer());
257 ret->gen_invert_this();
258 return ret;
259}
260
261int
262RefSCMatrix::nrow() const
263{
264 if (null()) return 0;
265 else return pointer()->nrow();
266}
267
268int
269RefSCMatrix::ncol() const
270{
271 if (null()) return 0;
272 else return pointer()->ncol();
273}
274
275RefSCDimension
276RefSCMatrix::rowdim() const
277{
278 if (null()) return 0;
279 else return pointer()->rowdim();
280}
281
282RefSCDimension
283RefSCMatrix::coldim() const
284{
285 if (null()) return 0;
286 else return pointer()->coldim();
287}
288
289Ref<SCMatrixKit>
290RefSCMatrix::kit() const
291{
292 if (null()) return 0;
293 else return pointer()->kit();
294}
295
296SCMatrixdouble
297RefSCMatrix::operator()(int i,int j) const
298{
299 return SCMatrixdouble(pointer(),i,j);
300}
301
302RefSCMatrix
303RefSCMatrix::clone() const
304{
305 RefSCMatrix r = kit()->matrix(rowdim(),coldim());
306 return r;
307}
308
309RefSCMatrix
310RefSCMatrix::get_subblock(int br, int er, int bc, int ec)
311{
312 require_nonnull();
313
314 RefSCMatrix ret = pointer()->get_subblock(br,er,bc,ec);
315 return ret;
316}
317
318void
319RefSCMatrix::assign_subblock(const RefSCMatrix& sb,
320 int br, int er, int bc, int ec, int sbr, int sbc)
321{
322 require_nonnull();
323 sb.require_nonnull();
324 pointer()->assign_subblock(sb.pointer(),br,er,bc,ec,sbr,sbc);
325}
326
327void
328RefSCMatrix::accumulate_subblock(const RefSCMatrix& sb,
329 int br, int er, int bc, int ec,
330 int sbr, int sbc)
331{
332 require_nonnull();
333 sb.require_nonnull();
334 pointer()->accumulate_subblock(sb.pointer(),br,er,bc,ec,sbr,sbc);
335}
336
337RefSCVector
338RefSCMatrix::get_row(int i) const
339{
340 require_nonnull();
341
342 RefSCVector ret = pointer()->get_row(i);
343 return ret;
344}
345
346RefSCVector
347RefSCMatrix::get_column(int i) const
348{
349 require_nonnull();
350
351 RefSCVector ret = pointer()->get_column(i);
352 return ret;
353}
354
355void
356RefSCMatrix::assign_row(const RefSCVector& v, int i) const
357{
358 require_nonnull();
359 v.require_nonnull();
360 pointer()->assign_row(v.pointer(),i);
361}
362
363void
364RefSCMatrix::assign_column(const RefSCVector& v, int i) const
365{
366 require_nonnull();
367 v.require_nonnull();
368 pointer()->assign_column(v.pointer(),i);
369}
370
371void
372RefSCMatrix::accumulate_row(const RefSCVector& v, int i) const
373{
374 require_nonnull();
375 v.require_nonnull();
376 pointer()->accumulate_row(v.pointer(),i);
377}
378
379void
380RefSCMatrix::accumulate_column(const RefSCVector& v, int i) const
381{
382 require_nonnull();
383 v.require_nonnull();
384 pointer()->accumulate_column(v.pointer(),i);
385}
386
387void
388RefSCMatrix::accumulate_product(const RefSCMatrix&a,const RefSCMatrix&b) const
389{
390 require_nonnull();
391 pointer()->accumulate_product(a.pointer(),b.pointer());
392}
393
394RefSCMatrix
395RefSCMatrix::copy() const
396{
397 if (null()) return 0;
398 RefSCMatrix v = kit()->matrix(rowdim(),coldim());
399 v.assign(*this);
400 return v;
401}
402
403void
404RefSCMatrix::randomize() const
405{
406 require_nonnull();
407 pointer()->randomize();
408}
409
410void
411RefSCMatrix::assign(const RefSCMatrix&a) const
412{
413 require_nonnull();
414 pointer()->assign(a.pointer());
415}
416
417void
418RefSCMatrix::assign(const double*v) const
419{
420 require_nonnull();
421 pointer()->assign(v);
422}
423
424void
425RefSCMatrix::assign(const double**v) const
426{
427 require_nonnull();
428 pointer()->assign(v);
429}
430
431void
432RefSCMatrix::convert(double*v) const
433{
434 require_nonnull();
435 pointer()->convert(v);
436}
437
438void
439RefSCMatrix::convert(double**v) const
440{
441 require_nonnull();
442 pointer()->convert(v);
443}
444
445void
446RefSCMatrix::scale(double a) const
447{
448 require_nonnull();
449 pointer()->scale(a);
450}
451
452void
453RefSCMatrix::assign(double a) const
454{
455 require_nonnull();
456 pointer()->assign(a);
457}
458
459void
460RefSCMatrix::accumulate(const RefSCMatrix&a) const
461{
462 require_nonnull();
463 pointer()->accumulate(a.pointer());
464}
465
466void
467RefSCMatrix::accumulate(const RefSymmSCMatrix&a) const
468{
469 require_nonnull();
470 pointer()->accumulate(a.pointer());
471}
472
473void
474RefSCMatrix::accumulate(const RefDiagSCMatrix&a) const
475{
476 require_nonnull();
477 pointer()->accumulate(a.pointer());
478}
479
480void
481RefSCMatrix::element_op(const Ref<SCElementOp>&op) const
482{
483 if (nonnull()) pointer()->element_op(op);
484}
485
486void
487RefSCMatrix::element_op(const Ref<SCElementOp2>& op,
488 const RefSCMatrix& m) const
489{
490 if (nonnull()) pointer()->element_op(op,m.pointer());
491}
492
493void
494RefSCMatrix::element_op(const Ref<SCElementOp3>& op,
495 const RefSCMatrix& m,
496 const RefSCMatrix& n) const
497{
498 if (nonnull()) pointer()->element_op(op,m.pointer(),n.pointer());
499}
500
501void
502RefSCMatrix::print(ostream& out) const
503{
504 print(0,out);
505}
506
507void
508RefSCMatrix::print(const char*title,ostream&out, int precision) const
509{
510 if (nonnull()) {
511 pointer()->print(title,out,precision);
512 }
513 else {
514 if (title) out << endl << title << endl;
515 out << "null matrix" << endl;
516 }
517}
518
519RefSCMatrix
520RefSCMatrix::operator *(double a) const
521{
522 RefSCMatrix r(copy());
523 r.scale(a);
524 return r;
525}
526
527void
528RefSCMatrix::svd(const RefSCMatrix &U,
529 const RefDiagSCMatrix &sigma,
530 const RefSCMatrix &V)
531{
532 require_nonnull();
533 RefSCMatrix c = clone();
534 c->assign(pointer());
535 c->svd_this(U.pointer(), sigma.pointer(), V.pointer());
536}
537
538double
539RefSCMatrix::solve_lin(const RefSCVector& v) const
540{
541 require_nonnull();
542 RefSCMatrix c = clone();
543 c->assign(pointer());
544 return c->solve_this(v.pointer());
545}
546
547double
548RefSCMatrix::determ() const
549{
550 require_nonnull();
551 RefSCMatrix c = clone();
552 c->assign(pointer());
553 return c->determ_this();
554}
555
556double
557RefSCMatrix::trace() const
558{
559 require_nonnull();
560 return pointer()->trace();
561}
562
563RefSCMatrix
564operator *(double a, const RefSCMatrix& v)
565{
566 return v*a;
567}
568
569void
570RefSCMatrix::accumulate_outer_product(const RefSCVector& v1,
571 const RefSCVector&v2) const
572{
573 require_nonnull();
574 pointer()->accumulate_outer_product(v1.pointer(),v2.pointer());
575}
576
577void
578RefSCMatrix::save(StateOut&s)
579{
580 if (null()) s.put(0);
581 else {
582 s.put(1);
583 pointer()->save(s);
584 }
585}
586
587void
588RefSCMatrix::restore(StateIn&s)
589{
590 int have_matrix;
591 s.get(have_matrix);
592 if (have_matrix && nonnull()) {
593 pointer()->restore(s);
594 }
595 else if (have_matrix) {
596 ExEnv::errn() << "RefSCMatrix::restore: matrix not properly initialized" << endl;
597 abort();
598 }
599 else {
600 clear();
601 }
602}
603
604int
605RefSCMatrix::nblock() const
606{
607 BlockedSCMatrix *b = dynamic_cast<BlockedSCMatrix*>(pointer());
608 if (b) return b->nblocks();
609 return 1;
610}
611
612RefSCMatrix
613RefSCMatrix::block(int i) const
614{
615 BlockedSCMatrix *b = dynamic_cast<BlockedSCMatrix*>(pointer());
616 if (b) return b->block(i);
617 return *this;
618}
619
620///////////////////////////////////////////////////////////////////
621// RefSymmSCMatrix members
622
623RefSymmSCMatrix::RefSymmSCMatrix()
624{
625}
626
627RefSymmSCMatrix::RefSymmSCMatrix (const RefSymmSCMatrix & o):
628 Ref<SymmSCMatrix> (o)
629{
630}
631
632RefSymmSCMatrix::RefSymmSCMatrix (SymmSCMatrix * o):
633 Ref<SymmSCMatrix> (o)
634{
635}
636
637RefSymmSCMatrix::~RefSymmSCMatrix ()
638{
639}
640
641RefSymmSCMatrix&
642RefSymmSCMatrix::operator=(SymmSCMatrix* cr)
643{
644 Ref<SymmSCMatrix>::operator=(cr);
645 return *this;
646}
647
648RefSymmSCMatrix&
649RefSymmSCMatrix::operator=(const RefSymmSCMatrix & c)
650{
651 Ref<SymmSCMatrix>::operator=(c);
652 return *this;
653}
654
655RefSymmSCMatrix::RefSymmSCMatrix(const RefSCDimension&a,
656 const Ref<SCMatrixKit>&k)
657{
658 assign_pointer(k->symmmatrix(a));
659}
660
661void
662RefSymmSCMatrix::set_element(int i, int j, double a) const
663{
664 require_nonnull();
665 pointer()->set_element(i,j,a);
666}
667
668void
669RefSymmSCMatrix::accumulate_element(int i, int j, double a) const
670{
671 require_nonnull();
672 pointer()->accumulate_element(i,j,a);
673}
674
675double
676RefSymmSCMatrix::get_element(int i, int j) const
677{
678 require_nonnull();
679 return pointer()->get_element(i,j);
680}
681
682RefSCMatrix
683RefSymmSCMatrix::get_subblock(int br, int er, int bc, int ec)
684{
685 require_nonnull();
686
687 RefSCMatrix ret = pointer()->get_subblock(br,er,bc,ec);
688 return ret;
689}
690
691RefSymmSCMatrix
692RefSymmSCMatrix::get_subblock(int br, int er)
693{
694 require_nonnull();
695
696 RefSymmSCMatrix ret = pointer()->get_subblock(br,er);
697 return ret;
698}
699
700void
701RefSymmSCMatrix::assign_subblock(const RefSCMatrix& sb,
702 int br, int er, int bc, int ec)
703{
704 require_nonnull();
705 sb.require_nonnull();
706 pointer()->assign_subblock(sb.pointer(),br,er,bc,ec);
707}
708
709void
710RefSymmSCMatrix::assign_subblock(const RefSymmSCMatrix& sb, int br, int er)
711{
712 require_nonnull();
713 sb.require_nonnull();
714 pointer()->assign_subblock(sb.pointer(),br,er);
715}
716
717void
718RefSymmSCMatrix::accumulate_subblock(const RefSCMatrix& sb,
719 int br, int er, int bc, int ec)
720{
721 require_nonnull();
722 sb.require_nonnull();
723 pointer()->accumulate_subblock(sb.pointer(),br,er,bc,ec);
724}
725
726void
727RefSymmSCMatrix::accumulate_subblock(const RefSymmSCMatrix& sb, int br, int er)
728{
729 require_nonnull();
730 sb.require_nonnull();
731 pointer()->accumulate_subblock(sb.pointer(),br,er);
732}
733
734RefSCVector
735RefSymmSCMatrix::get_row(int i)
736{
737 require_nonnull();
738
739 RefSCVector ret = pointer()->get_row(i);
740 return ret;
741}
742
743void
744RefSymmSCMatrix::assign_row(const RefSCVector& v, int i)
745{
746 require_nonnull();
747 v.require_nonnull();
748 pointer()->assign_row(v.pointer(),i);
749}
750
751void
752RefSymmSCMatrix::accumulate_row(const RefSCVector& v, int i)
753{
754 require_nonnull();
755 v.require_nonnull();
756 pointer()->accumulate_row(v.pointer(),i);
757}
758
759void
760RefSymmSCMatrix::accumulate_symmetric_product(const RefSCMatrix& a) const
761{
762 require_nonnull();
763 pointer()->accumulate_symmetric_product(a.pointer());
764}
765
766void
767RefSymmSCMatrix::accumulate_symmetric_sum(const RefSCMatrix& a) const
768{
769 require_nonnull();
770 pointer()->accumulate_symmetric_sum(a.pointer());
771}
772
773void
774RefSymmSCMatrix::accumulate_transform(const RefSCMatrix& a,
775 const RefSymmSCMatrix&b,
776 SCMatrix::Transform t) const
777{
778 require_nonnull();
779 pointer()->accumulate_transform(a.pointer(),b.pointer(),t);
780}
781
782void
783RefSymmSCMatrix::accumulate_transform(const RefSCMatrix& a,
784 const RefDiagSCMatrix&b,
785 SCMatrix::Transform t) const
786{
787 require_nonnull();
788 pointer()->accumulate_transform(a.pointer(),b.pointer(),t);
789}
790
791void
792RefSymmSCMatrix::accumulate_transform(const RefSymmSCMatrix& a,
793 const RefSymmSCMatrix&b) const
794{
795 require_nonnull();
796 pointer()->accumulate_transform(a.pointer(),b.pointer());
797}
798
799RefSymmSCMatrix
800RefSymmSCMatrix::operator+(const RefSymmSCMatrix&a) const
801{
802 require_nonnull();
803 a.require_nonnull();
804
805 RefSymmSCMatrix ret(dim(),kit());
806
807 ret->assign(pointer());
808 ret->accumulate(a.pointer());
809
810 return ret;
811}
812
813RefSymmSCMatrix
814RefSymmSCMatrix::operator-(const RefSymmSCMatrix&a) const
815{
816 require_nonnull();
817 a.require_nonnull();
818
819 RefSymmSCMatrix ret(dim(),kit());
820
821 ret->assign(a.pointer());
822 ret->scale(-1.0);
823 ret->accumulate(pointer());
824
825 return ret;
826}
827
828RefSymmSCMatrix
829RefSymmSCMatrix::i() const
830{
831 require_nonnull();
832
833 RefSymmSCMatrix ret;
834 ret = clone();
835 ret->assign(pointer());
836 ret->invert_this();
837 return ret;
838}
839
840RefSymmSCMatrix
841RefSymmSCMatrix::gi() const
842{
843 require_nonnull();
844
845 RefSymmSCMatrix ret;
846 ret = clone();
847 ret->assign(pointer());
848 ret->gen_invert_this();
849 return ret;
850}
851
852int
853RefSymmSCMatrix::n() const
854{
855 if (null()) return 0;
856 else return pointer()->dim()->n();
857}
858
859RefSCDimension
860RefSymmSCMatrix::dim() const
861{
862 if (null()) return 0;
863 else return pointer()->dim();
864}
865
866Ref<SCMatrixKit>
867RefSymmSCMatrix::kit() const
868{
869 if (null()) return 0;
870 else return pointer()->kit();
871}
872
873SymmSCMatrixdouble
874RefSymmSCMatrix::operator()(int i,int j) const
875{
876 return SymmSCMatrixdouble(pointer(),i,j);
877}
878
879RefSymmSCMatrix
880RefSymmSCMatrix::clone() const
881{
882 RefSymmSCMatrix r = kit()->symmmatrix(dim());
883 return r;
884}
885
886RefDiagSCMatrix
887RefSymmSCMatrix::eigvals() const
888{
889 if (null()) return 0;
890 RefDiagSCMatrix vals = kit()->diagmatrix(dim());
891 RefSCMatrix vecs = kit()->matrix(dim(),dim());
892 diagonalize(vals,vecs);
893 return vals;
894}
895
896RefSCMatrix
897RefSymmSCMatrix::eigvecs() const
898{
899 if (null()) return 0;
900 RefDiagSCMatrix vals = kit()->diagmatrix(dim());
901 RefSCMatrix vecs = kit()->matrix(dim(),dim());
902 diagonalize(vals,vecs);
903 return vecs;
904}
905
906double
907RefSymmSCMatrix::solve_lin(const RefSCVector& v) const
908{
909 require_nonnull();
910
911 RefSymmSCMatrix ret = clone();
912 ret->assign(pointer());
913 return ret->solve_this(v.pointer());
914}
915
916double
917RefSymmSCMatrix::determ() const
918{
919 require_nonnull();
920
921 RefSymmSCMatrix ret = clone();
922 ret->assign(pointer());
923 return ret->determ_this();
924}
925
926double
927RefSymmSCMatrix::trace() const
928{
929 require_nonnull();
930 return pointer()->trace();
931}
932
933void
934RefSymmSCMatrix::diagonalize(const RefDiagSCMatrix& vals,
935 const RefSCMatrix& vecs) const
936{
937 require_nonnull();
938 pointer()->diagonalize(vals.pointer(),vecs.pointer());
939}
940
941RefSymmSCMatrix
942RefSymmSCMatrix::copy() const
943{
944 if (null()) return 0;
945 RefSymmSCMatrix v = kit()->symmmatrix(dim());
946 v.assign(*this);
947 return v;
948}
949
950void
951RefSymmSCMatrix::randomize() const
952{
953 require_nonnull();
954 pointer()->randomize();
955}
956
957void
958RefSymmSCMatrix::assign(const RefSymmSCMatrix&a) const
959{
960 require_nonnull();
961 pointer()->assign(a.pointer());
962}
963
964void
965RefSymmSCMatrix::assign(const double*v) const
966{
967 require_nonnull();
968 pointer()->assign(v);
969}
970
971void
972RefSymmSCMatrix::assign(const double**v) const
973{
974 require_nonnull();
975 pointer()->assign(v);
976}
977
978void
979RefSymmSCMatrix::convert(double*v) const
980{
981 require_nonnull();
982 pointer()->convert(v);
983}
984
985void
986RefSymmSCMatrix::convert(double**v) const
987{
988 require_nonnull();
989 pointer()->convert(v);
990}
991
992void
993RefSymmSCMatrix::scale(double a) const
994{
995 require_nonnull();
996 pointer()->scale(a);
997}
998
999void
1000RefSymmSCMatrix::assign(double a) const
1001{
1002 require_nonnull();
1003 pointer()->assign(a);
1004}
1005
1006void
1007RefSymmSCMatrix::accumulate(const RefSymmSCMatrix&a) const
1008{
1009 require_nonnull();
1010 pointer()->accumulate(a.pointer());
1011}
1012
1013void
1014RefSymmSCMatrix::element_op(const Ref<SCElementOp>&op) const
1015{
1016 if (nonnull()) pointer()->element_op(op);
1017}
1018
1019void
1020RefSymmSCMatrix::element_op(const Ref<SCElementOp2>&op,
1021 const RefSymmSCMatrix&m) const
1022{
1023 if (nonnull()) pointer()->element_op(op,m.pointer());
1024}
1025
1026void
1027RefSymmSCMatrix::element_op(const Ref<SCElementOp3>&op,
1028 const RefSymmSCMatrix&m,
1029 const RefSymmSCMatrix&n) const
1030{
1031 if (nonnull()) pointer()->element_op(op,m.pointer(),n.pointer());
1032}
1033
1034void
1035RefSymmSCMatrix::print(ostream& out) const
1036{
1037 print(0,out);
1038}
1039
1040void
1041RefSymmSCMatrix::print(const char*title,ostream&out, int precision) const
1042{
1043 if (nonnull()) {
1044 pointer()->print(title,out,precision);
1045 }
1046 else {
1047 if (title) out << endl << title << endl;
1048 out << "null matrix" << endl;
1049 }
1050}
1051
1052RefSCMatrix
1053RefSymmSCMatrix::operator*(const RefSCMatrix&a) const
1054{
1055 require_nonnull();
1056 a.require_nonnull();
1057
1058 RefSCMatrix r = kit()->matrix(dim(),a->coldim());
1059 r->assign(0.0);
1060 r->accumulate_product(pointer(),a.pointer());
1061 return r;
1062}
1063
1064RefSCMatrix
1065RefSymmSCMatrix::operator*(const RefSymmSCMatrix&a) const
1066{
1067 require_nonnull();
1068 a.require_nonnull();
1069
1070 RefSCMatrix r = kit()->matrix(dim(),a->dim());
1071 r->assign(0.0);
1072 r->accumulate_product(pointer(),a.pointer());
1073 return r;
1074}
1075
1076RefSCVector
1077RefSymmSCMatrix::operator*(const RefSCVector&a) const
1078{
1079 require_nonnull();
1080 a.require_nonnull();
1081
1082 RefSCVector r = kit()->vector(dim());
1083 r->assign(0.0);
1084 r->accumulate_product(pointer(),a.pointer());
1085 return r;
1086}
1087
1088RefSymmSCMatrix
1089RefSymmSCMatrix::operator *(double a) const
1090{
1091 RefSymmSCMatrix r(copy());
1092 r.scale(a);
1093 return r;
1094}
1095
1096RefSymmSCMatrix
1097operator *(double a, const RefSymmSCMatrix& v)
1098{
1099 return v*a;
1100}
1101
1102void
1103RefSymmSCMatrix::accumulate_symmetric_outer_product(const RefSCVector&v) const
1104{
1105 require_nonnull();
1106 pointer()->accumulate_symmetric_outer_product(v.pointer());
1107}
1108
1109double
1110RefSymmSCMatrix::scalar_product(const RefSCVector&v) const
1111{
1112 if (null()) return 0.0;
1113 return pointer()->scalar_product(v.pointer());
1114}
1115
1116void
1117RefSymmSCMatrix::save(StateOut&s)
1118{
1119 if (null()) s.put(0);
1120 else {
1121 s.put(1);
1122 pointer()->save(s);
1123 }
1124}
1125
1126void
1127RefSymmSCMatrix::restore(StateIn&s)
1128{
1129 int have_matrix;
1130 s.get(have_matrix);
1131 if (have_matrix && nonnull()) {
1132 pointer()->restore(s);
1133 }
1134 else if (have_matrix) {
1135 ExEnv::errn() << "RefSymmSCMatrix::restore: "
1136 << "matrix not properly initialized" << endl;
1137 abort();
1138 }
1139 else {
1140 clear();
1141 }
1142}
1143
1144int
1145RefSymmSCMatrix::nblock() const
1146{
1147 BlockedSymmSCMatrix *b = dynamic_cast<BlockedSymmSCMatrix*>(pointer());
1148 if (b) return b->nblocks();
1149 return 1;
1150}
1151
1152RefSymmSCMatrix
1153RefSymmSCMatrix::block(int i) const
1154{
1155 BlockedSymmSCMatrix *b = dynamic_cast<BlockedSymmSCMatrix*>(pointer());
1156 if (b) return b->block(i);
1157 return *this;
1158}
1159
1160
1161///////////////////////////////////////////////////////////////////
1162// RefDiagSCMatrix members
1163
1164RefDiagSCMatrix::RefDiagSCMatrix()
1165{
1166}
1167
1168RefDiagSCMatrix::RefDiagSCMatrix (const RefDiagSCMatrix & o):
1169 Ref<DiagSCMatrix> (o)
1170{
1171}
1172
1173RefDiagSCMatrix::RefDiagSCMatrix (DiagSCMatrix * o):
1174 Ref<DiagSCMatrix> (o)
1175{
1176}
1177
1178RefDiagSCMatrix::~RefDiagSCMatrix ()
1179{
1180}
1181
1182RefDiagSCMatrix&
1183RefDiagSCMatrix::operator=(DiagSCMatrix* cr)
1184{
1185 Ref<DiagSCMatrix>::operator=(cr);
1186 return *this;
1187}
1188
1189RefDiagSCMatrix&
1190RefDiagSCMatrix::operator=(const RefDiagSCMatrix & c)
1191{
1192 Ref<DiagSCMatrix>::operator=(c);
1193 return *this;
1194}
1195
1196RefDiagSCMatrix::RefDiagSCMatrix(const RefSCDimension&a,
1197 const Ref<SCMatrixKit>&k)
1198{
1199 a.require_nonnull();
1200 assign_pointer(k->diagmatrix(a));
1201}
1202
1203void
1204RefDiagSCMatrix::set_element(int i, double a) const
1205{
1206 require_nonnull();
1207 pointer()->set_element(i,a);
1208}
1209
1210void
1211RefDiagSCMatrix::accumulate_element(int i, double a) const
1212{
1213 require_nonnull();
1214 pointer()->accumulate_element(i,a);
1215}
1216
1217double
1218RefDiagSCMatrix::get_element(int i) const
1219{
1220 require_nonnull();
1221 return pointer()->get_element(i);
1222}
1223
1224RefSCMatrix
1225RefDiagSCMatrix::operator*(const RefSCMatrix&a) const
1226{
1227 require_nonnull();
1228 a.require_nonnull();
1229
1230 RefSCMatrix r = kit()->matrix(dim(),a->coldim());
1231 r->assign(0.0);
1232 r->accumulate_product(pointer(),a.pointer());
1233 return r;
1234}
1235
1236RefDiagSCMatrix
1237RefDiagSCMatrix::operator+(const RefDiagSCMatrix&a) const
1238{
1239 require_nonnull();
1240 a.require_nonnull();
1241
1242 RefDiagSCMatrix ret(dim(),kit());
1243
1244 ret->assign(pointer());
1245 ret->accumulate(a.pointer());
1246
1247 return ret;
1248}
1249
1250RefDiagSCMatrix
1251RefDiagSCMatrix::operator-(const RefDiagSCMatrix&a) const
1252{
1253 require_nonnull();
1254 a.require_nonnull();
1255
1256 RefDiagSCMatrix ret(dim(),kit());
1257
1258 ret->assign(a.pointer());
1259 ret->scale(-1.0);
1260 ret->accumulate(pointer());
1261
1262 return ret;
1263}
1264
1265RefDiagSCMatrix
1266RefDiagSCMatrix::i() const
1267{
1268 require_nonnull();
1269
1270 RefDiagSCMatrix ret;
1271 ret = clone();
1272 ret->assign(pointer());
1273 ret->invert_this();
1274 return ret;
1275}
1276
1277RefDiagSCMatrix
1278RefDiagSCMatrix::gi() const
1279{
1280 require_nonnull();
1281
1282 RefDiagSCMatrix ret;
1283 ret = clone();
1284 ret->assign(pointer());
1285 ret->gen_invert_this();
1286 return ret;
1287}
1288
1289int
1290RefDiagSCMatrix::n() const
1291{
1292 if (null()) return 0;
1293 else return pointer()->dim()->n();
1294}
1295
1296RefSCDimension
1297RefDiagSCMatrix::dim() const
1298{
1299 if (null()) return 0;
1300 else return pointer()->dim();
1301}
1302
1303Ref<SCMatrixKit>
1304RefDiagSCMatrix::kit() const
1305{
1306 if (null()) return 0;
1307 else return pointer()->kit();
1308}
1309
1310DiagSCMatrixdouble
1311RefDiagSCMatrix::operator()(int i) const
1312{
1313 return DiagSCMatrixdouble(pointer(),i,i);
1314}
1315
1316RefDiagSCMatrix
1317RefDiagSCMatrix::clone() const
1318{
1319 RefDiagSCMatrix r = kit()->diagmatrix(dim());
1320 return r;
1321}
1322
1323RefDiagSCMatrix
1324RefDiagSCMatrix::copy() const
1325{
1326 if (null()) return 0;
1327 RefDiagSCMatrix v = kit()->diagmatrix(dim());
1328 v.assign(*this);
1329 return v;
1330}
1331
1332void
1333RefDiagSCMatrix::randomize() const
1334{
1335 require_nonnull();
1336 pointer()->randomize();
1337}
1338
1339void
1340RefDiagSCMatrix::assign(const RefDiagSCMatrix&a) const
1341{
1342 require_nonnull();
1343 pointer()->assign(a.pointer());
1344}
1345
1346void
1347RefDiagSCMatrix::assign(const double*v) const
1348{
1349 require_nonnull();
1350 pointer()->assign(v);
1351}
1352
1353void
1354RefDiagSCMatrix::convert(double*v) const
1355{
1356 require_nonnull();
1357 pointer()->convert(v);
1358}
1359
1360void
1361RefDiagSCMatrix::scale(double a) const
1362{
1363 require_nonnull();
1364 pointer()->scale(a);
1365}
1366
1367void
1368RefDiagSCMatrix::assign(double a) const
1369{
1370 require_nonnull();
1371 pointer()->assign(a);
1372}
1373
1374void
1375RefDiagSCMatrix::accumulate(const RefDiagSCMatrix&a) const
1376{
1377 require_nonnull();
1378 pointer()->accumulate(a.pointer());
1379}
1380
1381void
1382RefDiagSCMatrix::element_op(const Ref<SCElementOp>&op) const
1383{
1384 if (nonnull()) pointer()->element_op(op);
1385}
1386
1387void
1388RefDiagSCMatrix::element_op(const Ref<SCElementOp2>&op,
1389 const RefDiagSCMatrix&m) const
1390{
1391 if (nonnull()) pointer()->element_op(op,m.pointer());
1392}
1393
1394void
1395RefDiagSCMatrix::element_op(const Ref<SCElementOp3>&op,
1396 const RefDiagSCMatrix&m,
1397 const RefDiagSCMatrix&n) const
1398{
1399 if (nonnull()) pointer()->element_op(op,m.pointer(),n.pointer());
1400}
1401
1402double
1403RefDiagSCMatrix::determ() const
1404{
1405 return pointer()->determ_this();
1406}
1407
1408double
1409RefDiagSCMatrix::trace() const
1410{
1411 return pointer()->trace();
1412}
1413
1414void
1415RefDiagSCMatrix::print(ostream& out) const
1416{
1417 print(0,out);
1418}
1419
1420void
1421RefDiagSCMatrix::print(const char*title,ostream&out, int precision) const
1422{
1423 if (nonnull()) {
1424 pointer()->print(title,out,precision);
1425 }
1426 else {
1427 if (title) out << endl << title << endl;
1428 out << "null matrix" << endl;
1429 }
1430}
1431
1432RefDiagSCMatrix
1433RefDiagSCMatrix::operator *(double a) const
1434{
1435 RefDiagSCMatrix r(copy());
1436 r.scale(a);
1437 return r;
1438}
1439
1440RefDiagSCMatrix
1441operator *(double a, const RefDiagSCMatrix& v)
1442{
1443 return v*a;
1444}
1445
1446void
1447RefDiagSCMatrix::save(StateOut&s)
1448{
1449 if (null()) s.put(0);
1450 else {
1451 s.put(1);
1452 pointer()->save(s);
1453 }
1454}
1455
1456void
1457RefDiagSCMatrix::restore(StateIn&s)
1458{
1459 int have_matrix;
1460 s.get(have_matrix);
1461 if (have_matrix && nonnull()) {
1462 pointer()->restore(s);
1463 }
1464 else if (have_matrix) {
1465 ExEnv::errn() << "RefDiagSCMatrix::restore: "
1466 << "matrix not properly initialized" << endl;
1467 abort();
1468 }
1469 else {
1470 clear();
1471 }
1472}
1473
1474int
1475RefDiagSCMatrix::nblock() const
1476{
1477 BlockedDiagSCMatrix *b = dynamic_cast<BlockedDiagSCMatrix*>(pointer());
1478 if (b) return b->nblocks();
1479 return 1;
1480}
1481
1482RefDiagSCMatrix
1483RefDiagSCMatrix::block(int i) const
1484{
1485 BlockedDiagSCMatrix *b = dynamic_cast<BlockedDiagSCMatrix*>(pointer());
1486 if (b) return b->block(i);
1487 return *this;
1488}
1489
1490///////////////////////////////////////////////////////////////////
1491// RefSCVector members
1492
1493RefSCVector::RefSCVector()
1494{
1495}
1496
1497RefSCVector::RefSCVector (const RefSCVector & o):
1498 Ref<SCVector> (o)
1499{
1500}
1501
1502RefSCVector::RefSCVector (SCVector * o):
1503 Ref<SCVector> (o)
1504{
1505}
1506
1507RefSCVector::~RefSCVector ()
1508{
1509}
1510
1511RefSCVector&
1512RefSCVector::operator=(SCVector* cr)
1513{
1514 Ref<SCVector>::operator=(cr);
1515 return *this;
1516}
1517
1518RefSCVector&
1519RefSCVector::operator=(const RefSCVector & c)
1520{
1521 Ref<SCVector>::operator=(c);
1522 return *this;
1523}
1524
1525RefSCVector::RefSCVector(const RefSCDimension&a,
1526 const Ref<SCMatrixKit>&k)
1527{
1528 a.require_nonnull();
1529 assign_pointer(k->vector(a));
1530}
1531
1532void
1533RefSCVector::set_element(int i, double a) const
1534{
1535 require_nonnull();
1536 pointer()->set_element(i,a);
1537}
1538
1539void
1540RefSCVector::accumulate_element(int i, double a) const
1541{
1542 require_nonnull();
1543 pointer()->accumulate_element(i,a);
1544}
1545
1546double
1547RefSCVector::get_element(int i) const
1548{
1549 require_nonnull();
1550 return pointer()->get_element(i);
1551}
1552
1553RefSCVector
1554RefSCVector::operator+(const RefSCVector&a) const
1555{
1556 require_nonnull();
1557 a.require_nonnull();
1558
1559 RefSCVector ret(dim(),kit());
1560
1561 ret->assign(pointer());
1562 ret->accumulate(a.pointer());
1563
1564 return ret;
1565}
1566
1567RefSCVector
1568RefSCVector::operator-(const RefSCVector&a) const
1569{
1570 require_nonnull();
1571 a.require_nonnull();
1572
1573 RefSCVector ret(dim(),kit());
1574
1575 ret->assign(a.pointer());
1576 ret->scale(-1.0);
1577 ret->accumulate(pointer());
1578
1579 return ret;
1580}
1581
1582int
1583RefSCVector::n() const
1584{
1585 if (null()) return 0;
1586 else return pointer()->dim()->n();
1587}
1588
1589RefSCDimension
1590RefSCVector::dim() const
1591{
1592 if (null()) return 0;
1593 else return pointer()->dim();
1594}
1595
1596Ref<SCMatrixKit>
1597RefSCVector::kit() const
1598{
1599 if (null()) return 0;
1600 else return pointer()->kit();
1601}
1602
1603SCVectordouble
1604RefSCVector::operator()(int i) const
1605{
1606 return SCVectordouble(pointer(),i);
1607}
1608
1609SCVectordouble
1610RefSCVector::operator[](int i) const
1611{
1612 return SCVectordouble(pointer(),i);
1613}
1614
1615RefSCVector
1616RefSCVector::clone() const
1617{
1618 RefSCVector r = kit()->vector(dim());
1619 return r;
1620}
1621
1622RefSCVector
1623RefSCVector::copy() const
1624{
1625 if (null()) return 0;
1626 RefSCVector v = kit()->vector(dim());
1627 v.assign(*this);
1628 return v;
1629}
1630
1631double
1632RefSCVector::dot(const RefSCVector&a) const
1633{
1634 require_nonnull();
1635 return pointer()->scalar_product(a.pointer());
1636}
1637
1638double
1639RefSCVector::scalar_product(const RefSCVector&a) const
1640{
1641 require_nonnull();
1642 return pointer()->scalar_product(a.pointer());
1643}
1644
1645void
1646RefSCVector::randomize() const
1647{
1648 require_nonnull();
1649 pointer()->randomize();
1650}
1651
1652void
1653RefSCVector::assign(const RefSCVector&a) const
1654{
1655 require_nonnull();
1656 pointer()->assign(a.pointer());
1657}
1658
1659void
1660RefSCVector::assign(const double*v) const
1661{
1662 require_nonnull();
1663 pointer()->assign(v);
1664}
1665
1666void
1667RefSCVector::convert(double*v) const
1668{
1669 require_nonnull();
1670 pointer()->convert(v);
1671}
1672
1673void
1674RefSCVector::scale(double a) const
1675{
1676 require_nonnull();
1677 pointer()->scale(a);
1678}
1679
1680void
1681RefSCVector::assign(double a) const
1682{
1683 require_nonnull();
1684 pointer()->assign(a);
1685}
1686
1687void
1688RefSCVector::accumulate(const RefSCVector&a) const
1689{
1690 require_nonnull();
1691 pointer()->accumulate(a.pointer());
1692}
1693
1694void
1695RefSCVector::accumulate_product(const RefSymmSCMatrix&a, const RefSCVector&b)
1696{
1697 require_nonnull();
1698 pointer()->accumulate_product(a.pointer(), b.pointer());
1699}
1700
1701void
1702RefSCVector::accumulate_product(const RefSCMatrix&a, const RefSCVector&b)
1703{
1704 require_nonnull();
1705 pointer()->accumulate_product(a.pointer(), b.pointer());
1706}
1707
1708void
1709RefSCVector::element_op(const Ref<SCElementOp>&op) const
1710{
1711 if (nonnull()) pointer()->element_op(op);
1712}
1713
1714void
1715RefSCVector::element_op(const Ref<SCElementOp2>&op,
1716 const RefSCVector&v) const
1717{
1718 if (nonnull()) pointer()->element_op(op,v.pointer());
1719}
1720
1721void
1722RefSCVector::element_op(const Ref<SCElementOp3>&op,
1723 const RefSCVector&v,
1724 const RefSCVector&w) const
1725{
1726 if (nonnull()) pointer()->element_op(op,v.pointer(),w.pointer());
1727}
1728
1729void
1730RefSCVector::print(ostream& out) const
1731{
1732 print(0,out);
1733}
1734
1735void
1736RefSCVector::print(const char*title,ostream&out, int precision) const
1737{
1738 if (nonnull()) {
1739 pointer()->print(title,out,precision);
1740 }
1741 else {
1742 if (title) out << endl << title << endl;
1743 out << "null matrix" << endl;
1744 }
1745}
1746
1747RefSCVector
1748RefSCVector::operator *(double a) const
1749{
1750 RefSCVector r(copy());
1751 r.scale(a);
1752 return r;
1753}
1754
1755RefSCVector
1756operator *(double a, const RefSCVector& v)
1757{
1758 return v*a;
1759}
1760
1761void
1762RefSCVector::normalize() const
1763{
1764 require_nonnull();
1765 pointer()->normalize();
1766}
1767
1768RefSymmSCMatrix
1769RefSCVector::symmetric_outer_product() const
1770{
1771 RefSymmSCMatrix result(dim(),kit());
1772 result.assign(0.0);
1773 result.accumulate_symmetric_outer_product(pointer());
1774 return result;
1775}
1776
1777RefSCMatrix
1778RefSCVector::outer_product(const RefSCVector&v) const
1779{
1780 RefSCMatrix result(dim(),v.dim(),kit());
1781 result.assign(0.0);
1782 result.accumulate_outer_product(*this,v);
1783 return result;
1784}
1785
1786double
1787RefSCVector::maxabs() const
1788{
1789 if (null()) return 0.0;
1790 return pointer()->maxabs();
1791}
1792
1793
1794void
1795RefSCVector::save(StateOut&s)
1796{
1797 if (null()) s.put(0);
1798 else {
1799 s.put(1);
1800 pointer()->save(s);
1801 }
1802}
1803
1804void
1805RefSCVector::restore(StateIn&s)
1806{
1807 int have_matrix;
1808 s.get(have_matrix);
1809 if (have_matrix && nonnull()) {
1810 pointer()->restore(s);
1811 }
1812 else if (have_matrix) {
1813 ExEnv::errn() << "RefSCVector::restore: vector not properly initialized" << endl;
1814 abort();
1815 }
1816 else {
1817 clear();
1818 }
1819}
1820
1821/////////////////////////////////////////////////////////////////////////////
1822
1823}
1824
1825// Local Variables:
1826// mode: c++
1827// c-file-style: "CLJ"
1828// End:
Note: See TracBrowser for help on using the repository browser.