Mercurial > forge

comparison extra/linear-algebra/ov-re-tri.cc @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 276d676c91da
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 #include "ov-re-tri.h"
2
3
4 octave_tri::octave_tri(const Matrix &m, tri_type t):
5 octave_matrix(m), tri(t)
6 {
7 }
8
9 octave_tri::octave_tri (const octave_tri& T):
10 octave_matrix(T), tri(T.tri)
11 {
12 }
13
14 octave_tri::~octave_tri(void)
15 {
16 }
17
18 octave_value *octave_tri::clone(void)
19 {
20 return new octave_tri(*this);
21 }
22
23 static octave_value *
24 tri_numeric_conversion_function(const octave_value& a)
25 {
26 CAST_CONV_ARG (const octave_tri &);
27
28 return new octave_matrix (v.matrix_value());
29 }
30
31 type_conv_fcn
32 octave_tri::numeric_conversion_function (void) const
33 {
34 return tri_numeric_conversion_function;
35 }
36
37 octave_value * octave_tri::try_narrowing_conversion(void)
38 {
39 octave_value *retval = octave_matrix::try_narrowing_conversion();
40
41 if ( retval==0){
42 int nr = matrix.rows ();
43 int nc = matrix.cols ();
44
45 bool istri=true;
46
47 if(tri==Lower){
48 for(int r=0; r<nr; r++)
49 for(int c=r+1; c<nc; c++)
50 if(matrix(r,c)!=0.0){
51 istri=false;
52 break;
53 }
54 }
55 else if(tri==Upper){
56 for(int c=0; c<nc; c++)
57 for(int r=c+1; r<nr; r++)
58 if(matrix(r,c)!=0.0){
59 istri=false;
60 break;
61 }
62 }
63 else
64 error("Not a triangular Matrix");
65
66 if(!istri)
67 retval = new octave_matrix (matrix);
68 }
69
70 return retval;
71 }
72
73 void octave_tri::assign (const octave_value_list& idx, const Matrix& rhs)
74 {
75 octave_matrix::assign(idx, rhs);
76 return;
77 }
78
79 octave_value octave_tri::transpose(void) const
80 {
81 return new octave_tri(this->matrix_value().transpose(), tri_type(! bool(tri)));
82 }
83
84 void octave_tri::print (ostream& os, bool pr_as_read_syntax = false) const
85 {
86 octave_matrix::print(os, pr_as_read_syntax);
87 os << (tri == Upper ? "Upper" : "Lower") << " Triangular";
88 newline(os);
89 }
90
91 DEFUNOP (transpose, tri)
92 {
93 CAST_UNOP_ARG (const octave_tri&);
94 return v.transpose ();
95 }
96
97 DEFBINOP(ldiv, tri, matrix)
98 {
99 CAST_BINOP_ARGS (const octave_tri&, const octave_matrix&);
100 const Matrix X = v1.matrix_value();
101 const Matrix Y = v2.matrix_value();
102
103 if(X.cols()!=Y.rows()){
104 error("ldiv -- X.cols!=Y.rows");
105 return octave_value();
106 }
107
108 if(X.cols()!=X.rows()){
109 error("ldiv -- X not square matrix");
110 return octave_value();
111 }
112
113 Matrix A(X.rows(), Y.cols());
114
115 if (v1.tri_value() == octave_tri::Lower){
116 for(int c=0; c< A.cols(); c++)
117 for(int r=0; r<A.rows(); r++){
118 double sum=Y(r,c);
119 for(int i=0; i<r; i++)
120 sum=sum-X(r,i)*A(i,c);
121 A(r,c)=sum/X(r,r);
122 }
123 }
124 else if (v1.tri_value() == octave_tri::Upper){
125 for(int c=0; c< A.cols(); c++)
126 for(int r=A.rows()-1; r>=0; r--){
127 double sum=Y(r,c);
128 for(int i=r+1; i<A.rows(); i++)
129 sum=sum-X(r,i)*A(i,c);
130 A(r,c)=sum/X(r,r);
131 }
132 }
133 else{
134 error("Not a triangular matrix");
135 }
136
137 return A;
138 }
139
140 DEFASSIGNOP (assign, tri, matrix)
141 {
142 CAST_BINOP_ARGS (octave_tri &, octave_matrix&);
143 v1.assign(idx, v2.matrix_value());
144 return octave_value();
145 }
146
147 void install_tri_ops(void)
148 {
149 INSTALL_UNOP (op_transpose, octave_tri, transpose);
150 INSTALL_UNOP (op_hermitian, octave_tri, transpose);
151
152 INSTALL_BINOP (op_ldiv, octave_tri, octave_matrix, ldiv);
153 INSTALL_ASSIGNOP (op_asn_eq, octave_tri, octave_matrix, assign);
154 }
155
156
157 DEFINE_OCTAVE_ALLOCATOR (octave_tri);
158
159 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_tri, "tri");
160