1 //$$ submat.cpp                         submatrices
2 
3 // Copyright (C) 1991,2,3,4: R B Davies
4 
5 #include <ossim/matrix/include.h>
6 
7 #include <ossim/matrix/newmat.h>
8 #include <ossim/matrix/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) const23 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) const35 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) const45 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) const54 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) const64 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) const73 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 
operator <<(const BaseMatrix & bmx)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 
operator =(const BaseMatrix & bmx)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 
operator <<(const Real * r)161 void GetSubMatrix::operator<<(const Real* r)
162 {
163    REPORT
164    Tracer tr("SubMatrix(<<Real*)");
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 
operator <<(const int * r)178 void GetSubMatrix::operator<<(const int* r)
179 {
180    REPORT
181    Tracer tr("SubMatrix(<<int*)");
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 
operator =(Real r)195 void GetSubMatrix::operator=(Real r)
196 {
197    REPORT
198    Tracer tr("SubMatrix(=Real)");
199    SetUpLHS();
200    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
201                                   // do need LoadOnEntry
202    MatrixRowCol sub; int i = row_number;
203    while (i--)
204    {
205       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
206       sub.Copy(r); mr.Next();
207    }
208 }
209 
Inject(const GeneralMatrix & gmx)210 void GetSubMatrix::Inject(const GeneralMatrix& gmx)
211 {
212    REPORT
213    Tracer tr("SubMatrix(inject)");
214    SetUpLHS();
215    if (row_number != gmx.Nrows() || col_number != gmx.Ncols())
216       Throw(IncompatibleDimensionsException());
217    MatrixRow mrx((GeneralMatrix*)(&gmx), LoadOnEntry);
218    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
219                                   // do need LoadOnEntry
220    MatrixRowCol sub; int i = row_number;
221    while (i--)
222    {
223       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
224       sub.Inject(mrx); mr.Next(); mrx.Next();
225    }
226 }
227 
operator +=(const BaseMatrix & bmx)228 void GetSubMatrix::operator+=(const BaseMatrix& bmx)
229 {
230    REPORT
231    Tracer tr("SubMatrix(+=)"); GeneralMatrix* gmx = 0;
232    // MatrixConversionCheck mcc;         // Check for loss of info
233    Try
234    {
235       SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
236       if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
237          Throw(IncompatibleDimensionsException());
238       MatrixRow mrx(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.Check(mrx);                            // check for loss of info
246          sub.Add(mrx); mr.Next(); mrx.Next();
247       }
248       gmx->tDelete();
249    }
250 
251    CatchAll
252    {
253       if (gmx) gmx->tDelete();
254       ReThrow;
255    }
256 }
257 
operator -=(const BaseMatrix & bmx)258 void GetSubMatrix::operator-=(const BaseMatrix& bmx)
259 {
260    REPORT
261    Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
262    // MatrixConversionCheck mcc;         // Check for loss of info
263    Try
264    {
265       SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
266       if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
267          Throw(IncompatibleDimensionsException());
268       MatrixRow mrx(gmx, LoadOnEntry);
269       MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
270                                      // do need LoadOnEntry
271       MatrixRowCol sub; int i = row_number;
272       while (i--)
273       {
274          mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
275          sub.Check(mrx);                            // check for loss of info
276          sub.Sub(mrx); mr.Next(); mrx.Next();
277       }
278       gmx->tDelete();
279    }
280 
281    CatchAll
282    {
283       if (gmx) gmx->tDelete();
284       ReThrow;
285    }
286 }
287 
operator +=(Real r)288 void GetSubMatrix::operator+=(Real r)
289 {
290    REPORT
291    Tracer tr("SubMatrix(+= or -= Real)");
292    // MatrixConversionCheck mcc;         // Check for loss of info
293    Try
294    {
295       SetUpLHS();
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.Check();                               // check for loss of info
303          sub.Add(r); mr.Next();
304       }
305    }
306 
307    CatchAll
308    {
309       ReThrow;
310    }
311 }
312 
operator *=(Real r)313 void GetSubMatrix::operator*=(Real r)
314 {
315    REPORT
316    Tracer tr("SubMatrix(*= or /= Real)");
317    // MatrixConversionCheck mcc;         // Check for loss of info
318    Try
319    {
320       SetUpLHS();
321       MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
322                                      // do need LoadOnEntry
323       MatrixRowCol sub; int i = row_number;
324       while (i--)
325       {
326          mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
327          sub.Multiply(r); mr.Next();
328       }
329    }
330 
331    CatchAll
332    {
333       ReThrow;
334    }
335 }
336 
337 #ifdef use_namespace
338 }
339 #endif
340 
341