source: ntrip/trunk/BNC/newmat/submat.cpp@ 9012

Last change on this file since 9012 was 8901, checked in by stuerze, 5 years ago

upgrade to newmat11 library
File size: 11.7 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 if (gm->type().is_symmetric() &&
260 ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
261 Throw(ProgramException("Illegal operation on symmetric"));
262 MatrixRow mrx(gmx, LoadOnEntry);
263 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
264 // do need LoadOnEntry
265 MatrixRowCol sub; int i = row_number;
266 while (i--)
267 {
268 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
269 sub.Check(mrx); // check for loss of info
270 sub.Add(mrx); mr.Next(); mrx.Next();
271 }
272 gmx->tDelete();
273 }
274
275 CatchAll
276 {
277 if (gmx) gmx->tDelete();
278 ReThrow;
279 }
280}
281
282void GetSubMatrix::SP_eq(const BaseMatrix& bmx)
283{
284 REPORT
285 Tracer tr("SubMatrix(SP_eq)"); GeneralMatrix* gmx = 0;
286 // MatrixConversionCheck mcc; // Check for loss of info
287 Try
288 {
289 SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
290 if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
291 Throw(IncompatibleDimensionsException());
292 if (gm->type().is_symmetric() &&
293 ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
294 Throw(ProgramException("Illegal operation on symmetric"));
295 MatrixRow mrx(gmx, LoadOnEntry);
296 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
297 // do need LoadOnEntry
298 MatrixRowCol sub; int i = row_number;
299 while (i--)
300 {
301 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
302 sub.Multiply(mrx); mr.Next(); mrx.Next();
303 }
304 gmx->tDelete();
305 }
306
307 CatchAll
308 {
309 if (gmx) gmx->tDelete();
310 ReThrow;
311 }
312}
313
314void GetSubMatrix::operator-=(const BaseMatrix& bmx)
315{
316 REPORT
317 Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
318 // MatrixConversionCheck mcc; // Check for loss of info
319 Try
320 {
321 SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
322 if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
323 Throw(IncompatibleDimensionsException());
324 if (gm->type().is_symmetric() &&
325 ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
326 Throw(ProgramException("Illegal operation on symmetric"));
327 MatrixRow mrx(gmx, LoadOnEntry);
328 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
329 // do need LoadOnEntry
330 MatrixRowCol sub; int i = row_number;
331 while (i--)
332 {
333 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
334 sub.Check(mrx); // check for loss of info
335 sub.Sub(mrx); mr.Next(); mrx.Next();
336 }
337 gmx->tDelete();
338 }
339
340 CatchAll
341 {
342 if (gmx) gmx->tDelete();
343 ReThrow;
344 }
345}
346
347void GetSubMatrix::operator+=(Real r)
348{
349 REPORT
350 Tracer tr("SubMatrix(+= or -= Real)");
351 // MatrixConversionCheck mcc; // Check for loss of info
352 Try
353 {
354 SetUpLHS();
355 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
356 // do need LoadOnEntry
357 MatrixRowCol sub; int i = row_number;
358 while (i--)
359 {
360 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
361 sub.Check(); // check for loss of info
362 sub.Add(r); mr.Next();
363 }
364 }
365
366 CatchAll
367 {
368 ReThrow;
369 }
370}
371
372void GetSubMatrix::operator*=(Real r)
373{
374 REPORT
375 Tracer tr("SubMatrix(*= or /= Real)");
376 // MatrixConversionCheck mcc; // Check for loss of info
377 Try
378 {
379 SetUpLHS();
380 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
381 // do need LoadOnEntry
382 MatrixRowCol sub; int i = row_number;
383 while (i--)
384 {
385 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
386 sub.Multiply(r); mr.Next();
387 }
388 }
389
390 CatchAll
391 {
392 ReThrow;
393 }
394}
395
396#ifdef use_namespace
397}
398#endif
399
400///@}
401
Note: See TracBrowser for help on using the repository browser.