source: ntrip/branches/BNC_LM/newmat/submat.cpp@ 9183

Last change on this file since 9183 was 2013, checked in by mervart, 15 years ago

* empty log message *
File size: 10.0 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file submat.cpp
5/// Submatrix manipulation.
6
7// Copyright (C) 1991,2,3,4: R B Davies
8
9#include "include.h"
10
11#include "newmat.h"
12#include "newmatrc.h"
13
14#ifdef use_namespace
15namespace NEWMAT {
16#endif
17
18#ifdef DO_REPORT
19#define REPORT { static ExeCounter ExeCount(__LINE__,11); ++ExeCount; }
20#else
21#define REPORT {}
22#endif
23
24
25/****************************** submatrices *********************************/
26
27GetSubMatrix BaseMatrix::submatrix(int first_row, int last_row, int first_col,
28 int last_col) const
29{
30 REPORT
31 Tracer tr("submatrix");
32 int a = first_row - 1; int b = last_row - first_row + 1;
33 int c = first_col - 1; int d = last_col - first_col + 1;
34 if (a<0 || b<0 || c<0 || d<0) Throw(SubMatrixDimensionException());
35 // allow zero rows or columns
36 return GetSubMatrix(this, a, b, c, d, false);
37}
38
39GetSubMatrix BaseMatrix::sym_submatrix(int first_row, int last_row) const
40{
41 REPORT
42 Tracer tr("sym_submatrix");
43 int a = first_row - 1; int b = last_row - first_row + 1;
44 if (a<0 || b<0) Throw(SubMatrixDimensionException());
45 // allow zero rows or columns
46 return GetSubMatrix( this, a, b, a, b, true);
47}
48
49GetSubMatrix BaseMatrix::row(int first_row) const
50{
51 REPORT
52 Tracer tr("SubMatrix(row)");
53 int a = first_row - 1;
54 if (a<0) Throw(SubMatrixDimensionException());
55 return GetSubMatrix(this, a, 1, 0, -1, false);
56}
57
58GetSubMatrix BaseMatrix::rows(int first_row, int last_row) const
59{
60 REPORT
61 Tracer tr("SubMatrix(rows)");
62 int a = first_row - 1; int b = last_row - first_row + 1;
63 if (a<0 || b<0) Throw(SubMatrixDimensionException());
64 // allow zero rows or columns
65 return GetSubMatrix(this, a, b, 0, -1, false);
66}
67
68GetSubMatrix BaseMatrix::column(int first_col) const
69{
70 REPORT
71 Tracer tr("SubMatrix(column)");
72 int c = first_col - 1;
73 if (c<0) Throw(SubMatrixDimensionException());
74 return GetSubMatrix(this, 0, -1, c, 1, false);
75}
76
77GetSubMatrix BaseMatrix::columns(int first_col, int last_col) const
78{
79 REPORT
80 Tracer tr("SubMatrix(columns)");
81 int c = first_col - 1; int d = last_col - first_col + 1;
82 if (c<0 || d<0) Throw(SubMatrixDimensionException());
83 // allow zero rows or columns
84 return GetSubMatrix(this, 0, -1, c, d, false);
85}
86
87void GetSubMatrix::SetUpLHS()
88{
89 REPORT
90 Tracer tr("SubMatrix(LHS)");
91 const BaseMatrix* bm1 = bm;
92 GeneralMatrix* gm1 = ((BaseMatrix*&)bm)->Evaluate();
93 if ((BaseMatrix*)gm1!=bm1)
94 Throw(ProgramException("Invalid LHS"));
95 if (row_number < 0) row_number = gm1->Nrows();
96 if (col_number < 0) col_number = gm1->Ncols();
97 if (row_skip+row_number > gm1->Nrows()
98 || col_skip+col_number > gm1->Ncols())
99 Throw(SubMatrixDimensionException());
100}
101
102void GetSubMatrix::operator<<(const BaseMatrix& bmx)
103{
104 REPORT
105 Tracer tr("SubMatrix(<<)"); GeneralMatrix* gmx = 0;
106 Try
107 {
108 SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
109 if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
110 Throw(IncompatibleDimensionsException());
111 MatrixRow mrx(gmx, LoadOnEntry);
112 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
113 // do need LoadOnEntry
114 MatrixRowCol sub; int i = row_number;
115 while (i--)
116 {
117 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
118 sub.Copy(mrx); mr.Next(); mrx.Next();
119 }
120 gmx->tDelete();
121 }
122
123 CatchAll
124 {
125 if (gmx) gmx->tDelete();
126 ReThrow;
127 }
128}
129
130void GetSubMatrix::operator=(const BaseMatrix& bmx)
131{
132 REPORT
133 Tracer tr("SubMatrix(=)"); GeneralMatrix* gmx = 0;
134 // MatrixConversionCheck mcc; // Check for loss of info
135 Try
136 {
137 SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
138 if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
139 Throw(IncompatibleDimensionsException());
140 LoadAndStoreFlag lasf =
141 ( row_skip == col_skip
142 && gm->type().is_symmetric()
143 && gmx->type().is_symmetric() )
144 ? LoadOnEntry+DirectPart
145 : LoadOnEntry;
146 MatrixRow mrx(gmx, lasf);
147 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
148 // do need LoadOnEntry
149 MatrixRowCol sub; int i = row_number;
150 while (i--)
151 {
152 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
153 sub.CopyCheck(mrx); mr.Next(); mrx.Next();
154 }
155 gmx->tDelete();
156 }
157
158 CatchAll
159 {
160 if (gmx) gmx->tDelete();
161 ReThrow;
162 }
163}
164
165void GetSubMatrix::operator<<(const double* r)
166{
167 REPORT
168 Tracer tr("SubMatrix(<<double*)");
169 SetUpLHS();
170 if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
171 Throw(SubMatrixDimensionException());
172 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
173 // do need LoadOnEntry
174 MatrixRowCol sub; int i = row_number;
175 while (i--)
176 {
177 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
178 sub.Copy(r); mr.Next();
179 }
180}
181
182void GetSubMatrix::operator<<(const float* r)
183{
184 REPORT
185 Tracer tr("SubMatrix(<<float*)");
186 SetUpLHS();
187 if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
188 Throw(SubMatrixDimensionException());
189 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
190 // do need LoadOnEntry
191 MatrixRowCol sub; int i = row_number;
192 while (i--)
193 {
194 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
195 sub.Copy(r); mr.Next();
196 }
197}
198
199void GetSubMatrix::operator<<(const int* r)
200{
201 REPORT
202 Tracer tr("SubMatrix(<<int*)");
203 SetUpLHS();
204 if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
205 Throw(SubMatrixDimensionException());
206 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
207 // do need LoadOnEntry
208 MatrixRowCol sub; int i = row_number;
209 while (i--)
210 {
211 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
212 sub.Copy(r); mr.Next();
213 }
214}
215
216void GetSubMatrix::operator=(Real r)
217{
218 REPORT
219 Tracer tr("SubMatrix(=Real)");
220 SetUpLHS();
221 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
222 // do need LoadOnEntry
223 MatrixRowCol sub; int i = row_number;
224 while (i--)
225 {
226 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
227 sub.Copy(r); mr.Next();
228 }
229}
230
231void GetSubMatrix::inject(const GeneralMatrix& gmx)
232{
233 REPORT
234 Tracer tr("SubMatrix(inject)");
235 SetUpLHS();
236 if (row_number != gmx.Nrows() || col_number != gmx.Ncols())
237 Throw(IncompatibleDimensionsException());
238 MatrixRow mrx((GeneralMatrix*)(&gmx), LoadOnEntry);
239 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
240 // do need LoadOnEntry
241 MatrixRowCol sub; int i = row_number;
242 while (i--)
243 {
244 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
245 sub.Inject(mrx); mr.Next(); mrx.Next();
246 }
247}
248
249void GetSubMatrix::operator+=(const BaseMatrix& bmx)
250{
251 REPORT
252 Tracer tr("SubMatrix(+=)"); GeneralMatrix* gmx = 0;
253 // MatrixConversionCheck mcc; // Check for loss of info
254 Try
255 {
256 SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
257 if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
258 Throw(IncompatibleDimensionsException());
259 MatrixRow mrx(gmx, LoadOnEntry);
260 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
261 // do need LoadOnEntry
262 MatrixRowCol sub; int i = row_number;
263 while (i--)
264 {
265 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
266 sub.Check(mrx); // check for loss of info
267 sub.Add(mrx); mr.Next(); mrx.Next();
268 }
269 gmx->tDelete();
270 }
271
272 CatchAll
273 {
274 if (gmx) gmx->tDelete();
275 ReThrow;
276 }
277}
278
279void GetSubMatrix::operator-=(const BaseMatrix& bmx)
280{
281 REPORT
282 Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
283 // MatrixConversionCheck mcc; // Check for loss of info
284 Try
285 {
286 SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
287 if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
288 Throw(IncompatibleDimensionsException());
289 MatrixRow mrx(gmx, LoadOnEntry);
290 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
291 // do need LoadOnEntry
292 MatrixRowCol sub; int i = row_number;
293 while (i--)
294 {
295 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
296 sub.Check(mrx); // check for loss of info
297 sub.Sub(mrx); mr.Next(); mrx.Next();
298 }
299 gmx->tDelete();
300 }
301
302 CatchAll
303 {
304 if (gmx) gmx->tDelete();
305 ReThrow;
306 }
307}
308
309void GetSubMatrix::operator+=(Real r)
310{
311 REPORT
312 Tracer tr("SubMatrix(+= or -= Real)");
313 // MatrixConversionCheck mcc; // Check for loss of info
314 Try
315 {
316 SetUpLHS();
317 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
318 // do need LoadOnEntry
319 MatrixRowCol sub; int i = row_number;
320 while (i--)
321 {
322 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
323 sub.Check(); // check for loss of info
324 sub.Add(r); mr.Next();
325 }
326 }
327
328 CatchAll
329 {
330 ReThrow;
331 }
332}
333
334void GetSubMatrix::operator*=(Real r)
335{
336 REPORT
337 Tracer tr("SubMatrix(*= or /= Real)");
338 // MatrixConversionCheck mcc; // Check for loss of info
339 Try
340 {
341 SetUpLHS();
342 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
343 // do need LoadOnEntry
344 MatrixRowCol sub; int i = row_number;
345 while (i--)
346 {
347 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
348 sub.Multiply(r); mr.Next();
349 }
350 }
351
352 CatchAll
353 {
354 ReThrow;
355 }
356}
357
358#ifdef use_namespace
359}
360#endif
361
362///@}
363
Note: See TracBrowser for help on using the repository browser.