1 : //$$ newmat5.cpp Transpose, evaluate etc
2 :
3 : // Copyright (C) 1991,2,3,4: R B Davies
4 :
5 : //#define WANT_STREAM
6 :
7 : #include "include.h"
8 :
9 : #include "newmat.h"
10 : #include "newmatrc.h"
11 :
12 : #ifdef use_namespace
13 : namespace NEWMAT {
14 : #endif
15 :
16 :
17 : #ifdef DO_REPORT
18 : #define REPORT { static ExeCounter ExeCount(__LINE__,5); ++ExeCount; }
19 : #else
20 : #define REPORT {}
21 : #endif
22 :
23 :
24 : /************************ carry out operations ******************************/
25 :
26 :
27 0 : GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
28 : {
29 : GeneralMatrix* gm1;
30 :
31 0 : if (Compare(Type().t(),mt))
32 : {
33 : REPORT
34 0 : gm1 = mt.New(ncols_value,nrows_value,tm);
35 0 : for (int i=0; i<ncols_value; i++)
36 : {
37 0 : MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
38 0 : MatrixCol mc(this, mr.Data(), LoadOnEntry, i);
39 : }
40 : }
41 : else
42 : {
43 : REPORT
44 0 : gm1 = mt.New(ncols_value,nrows_value,tm);
45 0 : MatrixRow mr(this, LoadOnEntry);
46 0 : MatrixCol mc(gm1, StoreOnExit+DirectPart);
47 0 : int i = nrows_value;
48 0 : while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
49 : }
50 0 : tDelete(); gm1->ReleaseAndDelete(); return gm1;
51 : }
52 :
53 0 : GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
54 0 : { REPORT return Evaluate(mt); }
55 :
56 :
57 0 : GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
58 0 : { REPORT return Evaluate(mt); }
59 :
60 0 : GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
61 : {
62 : REPORT
63 0 : GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
64 0 : gmx->nrows_value = 1; gmx->ncols_value = gmx->storage = storage;
65 0 : return BorrowStore(gmx,mt);
66 : }
67 :
68 0 : GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
69 : {
70 : REPORT
71 0 : GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
72 0 : gmx->ncols_value = 1; gmx->nrows_value = gmx->storage = storage;
73 0 : return BorrowStore(gmx,mt);
74 : }
75 :
76 0 : GeneralMatrix* IdentityMatrix::Transpose(TransposedMatrix*, MatrixType mt)
77 0 : { REPORT return Evaluate(mt); }
78 :
79 0 : GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
80 : {
81 0 : if (Compare(this->Type(),mt)) { REPORT return this; }
82 : REPORT
83 0 : GeneralMatrix* gmx = mt.New(nrows_value,ncols_value,this);
84 0 : MatrixRow mr(this, LoadOnEntry);
85 0 : MatrixRow mrx(gmx, StoreOnExit+DirectPart);
86 0 : int i=nrows_value;
87 0 : while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
88 0 : tDelete(); gmx->ReleaseAndDelete(); return gmx;
89 : }
90 :
91 0 : GeneralMatrix* GenericMatrix::Evaluate(MatrixType mt)
92 0 : { REPORT return gm->Evaluate(mt); }
93 :
94 0 : GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
95 : {
96 0 : gm=((BaseMatrix*&)bm)->Evaluate();
97 0 : int nr=gm->Nrows(); int nc=gm->Ncols();
98 0 : Compare(gm->Type().AddEqualEl(),mt);
99 0 : if (!(mt==gm->Type()))
100 : {
101 : REPORT
102 0 : GeneralMatrix* gmx = mt.New(nr,nc,this);
103 0 : MatrixRow mr(gm, LoadOnEntry);
104 0 : MatrixRow mrx(gmx, StoreOnExit+DirectPart);
105 0 : while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
106 0 : gmx->ReleaseAndDelete(); gm->tDelete();
107 0 : return gmx;
108 : }
109 0 : else if (gm->reuse())
110 : {
111 0 : REPORT gm->Add(f);
112 0 : return gm;
113 : }
114 : else
115 : {
116 0 : REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
117 0 : gmy->ReleaseAndDelete(); gmy->Add(gm,f);
118 0 : return gmy;
119 : }
120 : }
121 :
122 0 : GeneralMatrix* NegShiftedMatrix::Evaluate(MatrixType mt)
123 : {
124 0 : gm=((BaseMatrix*&)bm)->Evaluate();
125 0 : int nr=gm->Nrows(); int nc=gm->Ncols();
126 0 : Compare(gm->Type().AddEqualEl(),mt);
127 0 : if (!(mt==gm->Type()))
128 : {
129 : REPORT
130 0 : GeneralMatrix* gmx = mt.New(nr,nc,this);
131 0 : MatrixRow mr(gm, LoadOnEntry);
132 0 : MatrixRow mrx(gmx, StoreOnExit+DirectPart);
133 0 : while (nr--) { mrx.NegAdd(mr,f); mrx.Next(); mr.Next(); }
134 0 : gmx->ReleaseAndDelete(); gm->tDelete();
135 0 : return gmx;
136 : }
137 0 : else if (gm->reuse())
138 : {
139 0 : REPORT gm->NegAdd(f);
140 0 : return gm;
141 : }
142 : else
143 : {
144 0 : REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
145 0 : gmy->ReleaseAndDelete(); gmy->NegAdd(gm,f);
146 0 : return gmy;
147 : }
148 : }
149 :
150 0 : GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
151 : {
152 0 : gm=((BaseMatrix*&)bm)->Evaluate();
153 0 : int nr=gm->Nrows(); int nc=gm->Ncols();
154 0 : if (Compare(gm->Type(),mt))
155 : {
156 0 : if (gm->reuse())
157 : {
158 0 : REPORT gm->Multiply(f);
159 0 : return gm;
160 : }
161 : else
162 : {
163 0 : REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
164 0 : gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
165 0 : return gmx;
166 : }
167 : }
168 : else
169 : {
170 : REPORT
171 0 : GeneralMatrix* gmx = mt.New(nr,nc,this);
172 0 : MatrixRow mr(gm, LoadOnEntry);
173 0 : MatrixRow mrx(gmx, StoreOnExit+DirectPart);
174 0 : while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
175 0 : gmx->ReleaseAndDelete(); gm->tDelete();
176 0 : return gmx;
177 : }
178 : }
179 :
180 0 : GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
181 : {
182 0 : gm=((BaseMatrix*&)bm)->Evaluate();
183 0 : int nr=gm->Nrows(); int nc=gm->Ncols();
184 0 : if (Compare(gm->Type(),mt))
185 : {
186 0 : if (gm->reuse())
187 : {
188 0 : REPORT gm->Negate();
189 0 : return gm;
190 : }
191 : else
192 : {
193 : REPORT
194 0 : GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
195 0 : gmx->ReleaseAndDelete(); gmx->Negate(gm);
196 0 : return gmx;
197 : }
198 : }
199 : else
200 : {
201 : REPORT
202 0 : GeneralMatrix* gmx = mt.New(nr,nc,this);
203 0 : MatrixRow mr(gm, LoadOnEntry);
204 0 : MatrixRow mrx(gmx, StoreOnExit+DirectPart);
205 0 : while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
206 0 : gmx->ReleaseAndDelete(); gm->tDelete();
207 0 : return gmx;
208 : }
209 : }
210 :
211 0 : GeneralMatrix* ReversedMatrix::Evaluate(MatrixType mt)
212 : {
213 0 : gm=((BaseMatrix*&)bm)->Evaluate(); GeneralMatrix* gmx;
214 :
215 0 : if ((gm->Type()).IsBand() && ! (gm->Type()).IsDiagonal())
216 : {
217 0 : gm->tDelete();
218 0 : Throw(NotDefinedException("Reverse", "band matrices"));
219 : }
220 :
221 0 : if (gm->reuse()) { REPORT gm->ReverseElements(); gmx = gm; }
222 : else
223 : {
224 : REPORT
225 0 : gmx = gm->Type().New(gm->Nrows(), gm->Ncols(), this);
226 0 : gmx->ReverseElements(gm); gmx->ReleaseAndDelete();
227 : }
228 0 : return gmx->Evaluate(mt); // target matrix is different type?
229 :
230 : }
231 :
232 0 : GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
233 : {
234 : REPORT
235 0 : gm=((BaseMatrix*&)bm)->Evaluate();
236 0 : Compare(gm->Type().t(),mt);
237 0 : GeneralMatrix* gmx=gm->Transpose(this, mt);
238 0 : return gmx;
239 : }
240 :
241 0 : GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
242 : {
243 0 : gm = ((BaseMatrix*&)bm)->Evaluate();
244 0 : GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
245 0 : gmx->nrows_value = 1; gmx->ncols_value = gmx->storage = gm->storage;
246 0 : return gm->BorrowStore(gmx,mt);
247 : }
248 :
249 0 : GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
250 : {
251 0 : gm = ((BaseMatrix*&)bm)->Evaluate();
252 0 : GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
253 0 : gmx->ncols_value = 1; gmx->nrows_value = gmx->storage = gm->storage;
254 0 : return gm->BorrowStore(gmx,mt);
255 : }
256 :
257 0 : GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
258 : {
259 0 : gm = ((BaseMatrix*&)bm)->Evaluate();
260 0 : GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
261 0 : gmx->nrows_value = gmx->ncols_value = gmx->storage = gm->storage;
262 0 : return gm->BorrowStore(gmx,mt);
263 : }
264 :
265 0 : GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
266 : {
267 0 : Tracer tr("MatedMatrix::Evaluate");
268 0 : gm = ((BaseMatrix*&)bm)->Evaluate();
269 0 : GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
270 0 : gmx->nrows_value = nr; gmx->ncols_value = nc; gmx->storage = gm->storage;
271 0 : if (nr*nc != gmx->storage)
272 0 : Throw(IncompatibleDimensionsException());
273 0 : return gm->BorrowStore(gmx,mt);
274 : }
275 :
276 0 : GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
277 : {
278 : REPORT
279 0 : Tracer tr("SubMatrix(evaluate)");
280 0 : gm = ((BaseMatrix*&)bm)->Evaluate();
281 0 : if (row_number < 0) row_number = gm->Nrows();
282 0 : if (col_number < 0) col_number = gm->Ncols();
283 0 : if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
284 : {
285 0 : gm->tDelete();
286 0 : Throw(SubMatrixDimensionException());
287 : }
288 0 : if (IsSym) Compare(gm->Type().ssub(), mt);
289 0 : else Compare(gm->Type().sub(), mt);
290 0 : GeneralMatrix* gmx = mt.New(row_number, col_number, this);
291 0 : int i = row_number;
292 0 : MatrixRow mr(gm, LoadOnEntry, row_skip);
293 0 : MatrixRow mrx(gmx, StoreOnExit+DirectPart);
294 0 : MatrixRowCol sub;
295 0 : while (i--)
296 : {
297 0 : mr.SubRowCol(sub, col_skip, col_number); // put values in sub
298 0 : mrx.Copy(sub); mrx.Next(); mr.Next();
299 : }
300 0 : gmx->ReleaseAndDelete(); gm->tDelete();
301 0 : return gmx;
302 : }
303 :
304 :
305 0 : GeneralMatrix* ReturnMatrix::Evaluate(MatrixType mt)
306 : {
307 0 : return gm->Evaluate(mt);
308 : }
309 :
310 :
311 0 : void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
312 : {
313 : REPORT
314 0 : Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
315 0 : while (i--)
316 0 : { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
317 0 : i = storage & 3; while (i--) *s++ = *s1++ + f;
318 0 : }
319 :
320 0 : void GeneralMatrix::Add(Real f)
321 : {
322 : REPORT
323 0 : Real* s=store; int i=(storage >> 2);
324 0 : while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
325 0 : i = storage & 3; while (i--) *s++ += f;
326 0 : }
327 :
328 0 : void GeneralMatrix::NegAdd(GeneralMatrix* gm1, Real f)
329 : {
330 : REPORT
331 0 : Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
332 0 : while (i--)
333 0 : { *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; }
334 0 : i = storage & 3; while (i--) *s++ = f - *s1++;
335 0 : }
336 :
337 0 : void GeneralMatrix::NegAdd(Real f)
338 : {
339 : REPORT
340 0 : Real* s=store; int i=(storage >> 2);
341 0 : while (i--)
342 : {
343 0 : *s = f - *s; s++; *s = f - *s; s++;
344 0 : *s = f - *s; s++; *s = f - *s; s++;
345 : }
346 0 : i = storage & 3; while (i--) { *s = f - *s; s++; }
347 0 : }
348 :
349 0 : void GeneralMatrix::Negate(GeneralMatrix* gm1)
350 : {
351 : // change sign of elements
352 : REPORT
353 0 : Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
354 0 : while (i--)
355 0 : { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
356 0 : i = storage & 3; while(i--) *s++ = -(*s1++);
357 0 : }
358 :
359 0 : void GeneralMatrix::Negate()
360 : {
361 : REPORT
362 0 : Real* s=store; int i=(storage >> 2);
363 0 : while (i--)
364 0 : { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
365 0 : i = storage & 3; while(i--) { *s = -(*s); s++; }
366 0 : }
367 :
368 0 : void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
369 : {
370 : REPORT
371 0 : Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
372 0 : while (i--)
373 0 : { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
374 0 : i = storage & 3; while (i--) *s++ = *s1++ * f;
375 0 : }
376 :
377 0 : void GeneralMatrix::Multiply(Real f)
378 : {
379 : REPORT
380 0 : Real* s=store; int i=(storage >> 2);
381 0 : while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
382 0 : i = storage & 3; while (i--) *s++ *= f;
383 0 : }
384 :
385 :
386 : /************************ MatrixInput routines ****************************/
387 :
388 : // int MatrixInput::n; // number values still to be read
389 : // Real* MatrixInput::r; // pointer to next location to be read to
390 :
391 0 : MatrixInput MatrixInput::operator<<(double f)
392 : {
393 : REPORT
394 0 : Tracer et("MatrixInput");
395 0 : if (n<=0) Throw(ProgramException("List of values too long"));
396 0 : *r = (Real)f; int n1 = n-1; n=0; // n=0 so we won't trigger exception
397 0 : return MatrixInput(n1, r+1);
398 : }
399 :
400 :
401 0 : MatrixInput GeneralMatrix::operator<<(double f)
402 : {
403 : REPORT
404 0 : Tracer et("MatrixInput");
405 0 : int n = Storage();
406 0 : if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
407 0 : Real* r; r = Store(); *r = (Real)f; n--;
408 0 : return MatrixInput(n, r+1);
409 : }
410 :
411 0 : MatrixInput GetSubMatrix::operator<<(double f)
412 : {
413 : REPORT
414 0 : Tracer et("MatrixInput (GetSubMatrix)");
415 0 : SetUpLHS();
416 0 : if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
417 : {
418 0 : Throw(ProgramException("MatrixInput requires complete rows"));
419 : }
420 0 : MatrixRow mr(gm, DirectPart, row_skip); // to pick up location and length
421 0 : int n = mr.Storage();
422 0 : if (n<=0)
423 : {
424 0 : Throw(ProgramException("Loading data to zero length row"));
425 : }
426 0 : Real* r; r = mr.Data(); *r = (Real)f; n--;
427 0 : if (+(mr.cw*HaveStore))
428 : {
429 0 : Throw(ProgramException("Fails with this matrix type"));
430 : }
431 0 : return MatrixInput(n, r+1);
432 : }
433 :
434 0 : MatrixInput MatrixInput::operator<<(float f)
435 : {
436 : REPORT
437 0 : Tracer et("MatrixInput");
438 0 : if (n<=0) Throw(ProgramException("List of values too long"));
439 0 : *r = (Real)f; int n1 = n-1; n=0; // n=0 so we won't trigger exception
440 0 : return MatrixInput(n1, r+1);
441 : }
442 :
443 :
444 0 : MatrixInput GeneralMatrix::operator<<(float f)
445 : {
446 : REPORT
447 0 : Tracer et("MatrixInput");
448 0 : int n = Storage();
449 0 : if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
450 0 : Real* r; r = Store(); *r = (Real)f; n--;
451 0 : return MatrixInput(n, r+1);
452 : }
453 :
454 0 : MatrixInput GetSubMatrix::operator<<(float f)
455 : {
456 : REPORT
457 0 : Tracer et("MatrixInput (GetSubMatrix)");
458 0 : SetUpLHS();
459 0 : if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
460 : {
461 0 : Throw(ProgramException("MatrixInput requires complete rows"));
462 : }
463 0 : MatrixRow mr(gm, DirectPart, row_skip); // to pick up location and length
464 0 : int n = mr.Storage();
465 0 : if (n<=0)
466 : {
467 0 : Throw(ProgramException("Loading data to zero length row"));
468 : }
469 0 : Real* r; r = mr.Data(); *r = (Real)f; n--;
470 0 : if (+(mr.cw*HaveStore))
471 : {
472 0 : Throw(ProgramException("Fails with this matrix type"));
473 : }
474 0 : return MatrixInput(n, r+1);
475 : }
476 0 : MatrixInput::~MatrixInput()
477 : {
478 : REPORT
479 0 : Tracer et("MatrixInput");
480 0 : if (n!=0) Throw(ProgramException("A list of values was too short"));
481 0 : }
482 :
483 0 : MatrixInput BandMatrix::operator<<(double)
484 : {
485 0 : Tracer et("MatrixInput");
486 0 : bool dummy = true;
487 0 : if (dummy) // get rid of warning message
488 0 : Throw(ProgramException("Cannot use list read with a BandMatrix"));
489 0 : return MatrixInput(0, 0);
490 : }
491 :
492 0 : MatrixInput BandMatrix::operator<<(float)
493 : {
494 0 : Tracer et("MatrixInput");
495 0 : bool dummy = true;
496 0 : if (dummy) // get rid of warning message
497 0 : Throw(ProgramException("Cannot use list read with a BandMatrix"));
498 0 : return MatrixInput(0, 0);
499 : }
500 :
501 0 : void BandMatrix::operator<<(const double*)
502 0 : { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
503 :
504 0 : void BandMatrix::operator<<(const float*)
505 0 : { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
506 :
507 0 : void BandMatrix::operator<<(const int*)
508 0 : { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
509 :
510 0 : void SymmetricBandMatrix::operator<<(const double*)
511 0 : { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
512 :
513 0 : void SymmetricBandMatrix::operator<<(const float*)
514 0 : { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
515 :
516 0 : void SymmetricBandMatrix::operator<<(const int*)
517 0 : { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
518 :
519 : // ************************* Reverse order of elements ***********************
520 :
521 0 : void GeneralMatrix::ReverseElements(GeneralMatrix* gm)
522 : {
523 : // reversing into a new matrix
524 : REPORT
525 0 : int n = Storage(); Real* rx = Store() + n; Real* x = gm->Store();
526 0 : while (n--) *(--rx) = *(x++);
527 0 : }
528 :
529 0 : void GeneralMatrix::ReverseElements()
530 : {
531 : // reversing in place
532 : REPORT
533 0 : int n = Storage(); Real* x = Store(); Real* rx = x + n;
534 0 : n /= 2;
535 0 : while (n--) { Real t = *(--rx); *rx = *x; *(x++) = t; }
536 0 : }
537 :
538 :
539 : #ifdef use_namespace
540 78 : }
541 39 : #endif
542 :
|
