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