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
15 namespace 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
submatrix(int first_row,int last_row,int first_col,int last_col) const27 GetSubMatrix 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
sym_submatrix(int first_row,int last_row) const39 GetSubMatrix 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
row(int first_row) const49 GetSubMatrix 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
rows(int first_row,int last_row) const58 GetSubMatrix 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
column(int first_col) const68 GetSubMatrix 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
columns(int first_col,int last_col) const77 GetSubMatrix 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
SetUpLHS()87 void 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
operator <<(const BaseMatrix & bmx)102 void 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
operator =(const BaseMatrix & bmx)130 void 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
operator <<(const double * r)165 void 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
operator <<(const float * r)182 void 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
operator <<(const int * r)199 void 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
operator =(Real r)216 void 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
inject(const GeneralMatrix & gmx)231 void 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
operator +=(const BaseMatrix & bmx)249 void 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
SP_eq(const BaseMatrix & bmx)282 void 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
operator -=(const BaseMatrix & bmx)314 void 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
operator +=(Real r)347 void 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
operator *=(Real r)372 void 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
402