source: ntrip/trunk/BNS/newmat/newmat3.cpp@ 6183

Last change on this file since 6183 was 810, checked in by mervart, 17 years ago

* empty log message *
File size: 24.4 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file newmat3.cpp
5/// Get and restore rows and columns.
6
7// Copyright (C) 1991,2,3,4: R B Davies
8
9//#define WANT_STREAM
10
11#include "include.h"
12
13#include "newmat.h"
14#include "newmatrc.h"
15
16#ifdef use_namespace
17namespace NEWMAT {
18#endif
19
20
21#ifdef DO_REPORT
22#define REPORT { static ExeCounter ExeCount(__LINE__,3); ++ExeCount; }
23#else
24#define REPORT {}
25#endif
26
27//#define MONITOR(what,storage,store)
28// { cout << what << " " << storage << " at " << (long)store << "\n"; }
29
30#define MONITOR(what,store,storage) {}
31
32
33// Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
34//
35// LoadOnEntry:
36// Load data into MatrixRow or Col dummy array under GetRow or GetCol
37// StoreOnExit:
38// Restore data to original matrix under RestoreRow or RestoreCol
39// DirectPart:
40// Load or restore only part directly stored; must be set with StoreOnExit
41// Still have decide how to handle this with symmetric
42// StoreHere:
43// used in columns only - store data at supplied storage address;
44// used for GetCol, NextCol & RestoreCol. No need to fill out zeros
45// HaveStore:
46// dummy array has been assigned (internal use only).
47
48// For symmetric matrices, treat columns as rows unless StoreHere is set;
49// then stick to columns as this will give better performance for doing
50// inverses
51
52// How components are used:
53
54// Use rows wherever possible in preference to columns
55
56// Columns without StoreHere are used in in-exact transpose, sum column,
57// multiply a column vector, and maybe in future access to column,
58// additional multiply functions, add transpose
59
60// Columns with StoreHere are used in exact transpose (not symmetric matrices
61// or vectors, load only)
62
63// Columns with MatrixColX (Store to full column) are used in inverse and solve
64
65// Functions required for each matrix class
66
67// GetRow(MatrixRowCol& mrc)
68// GetCol(MatrixRowCol& mrc)
69// GetCol(MatrixColX& mrc)
70// RestoreRow(MatrixRowCol& mrc)
71// RestoreCol(MatrixRowCol& mrc)
72// RestoreCol(MatrixColX& mrc)
73// NextRow(MatrixRowCol& mrc)
74// NextCol(MatrixRowCol& mrc)
75// NextCol(MatrixColX& mrc)
76
77// The Restore routines assume StoreOnExit has already been checked
78// Defaults for the Next routines are given below
79// Assume cannot have both !DirectPart && StoreHere for MatrixRowCol routines
80
81
82// Default NextRow and NextCol:
83// will work as a default but need to override NextRow for efficiency
84
85void GeneralMatrix::NextRow(MatrixRowCol& mrc)
86{
87 REPORT
88 if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
89 mrc.rowcol++;
90 if (mrc.rowcol<nrows_val) { REPORT this->GetRow(mrc); }
91 else { REPORT mrc.cw -= StoreOnExit; }
92}
93
94void GeneralMatrix::NextCol(MatrixRowCol& mrc)
95{
96 REPORT // 423
97 if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
98 mrc.rowcol++;
99 if (mrc.rowcol<ncols_val) { REPORT this->GetCol(mrc); }
100 else { REPORT mrc.cw -= StoreOnExit; }
101}
102
103void GeneralMatrix::NextCol(MatrixColX& mrc)
104{
105 REPORT // 423
106 if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
107 mrc.rowcol++;
108 if (mrc.rowcol<ncols_val) { REPORT this->GetCol(mrc); }
109 else { REPORT mrc.cw -= StoreOnExit; }
110}
111
112
113// routines for matrix
114
115void Matrix::GetRow(MatrixRowCol& mrc)
116{
117 REPORT
118 mrc.skip=0; mrc.storage=mrc.length=ncols_val;
119 mrc.data=store+mrc.rowcol*ncols_val;
120}
121
122
123void Matrix::GetCol(MatrixRowCol& mrc)
124{
125 REPORT
126 mrc.skip=0; mrc.storage=mrc.length=nrows_val;
127 if ( ncols_val==1 && !(mrc.cw*StoreHere) ) // ColumnVector
128 { REPORT mrc.data=store; }
129 else
130 {
131 Real* ColCopy;
132 if ( !(mrc.cw*(HaveStore+StoreHere)) )
133 {
134 REPORT
135 ColCopy = new Real [nrows_val]; MatrixErrorNoSpace(ColCopy);
136 MONITOR_REAL_NEW("Make (MatGetCol)",nrows_val,ColCopy)
137 mrc.data = ColCopy; mrc.cw += HaveStore;
138 }
139 else { REPORT ColCopy = mrc.data; }
140 if (+(mrc.cw*LoadOnEntry))
141 {
142 REPORT
143 Real* Mstore = store+mrc.rowcol; int i=nrows_val;
144 //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
145 if (i) for (;;)
146 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
147 }
148 }
149}
150
151void Matrix::GetCol(MatrixColX& mrc)
152{
153 REPORT
154 mrc.skip=0; mrc.storage=nrows_val; mrc.length=nrows_val;
155 if (+(mrc.cw*LoadOnEntry))
156 {
157 REPORT Real* ColCopy = mrc.data;
158 Real* Mstore = store+mrc.rowcol; int i=nrows_val;
159 //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
160 if (i) for (;;)
161 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
162 }
163}
164
165void Matrix::RestoreCol(MatrixRowCol& mrc)
166{
167 // always check StoreOnExit before calling RestoreCol
168 REPORT // 429
169 if (+(mrc.cw*HaveStore))
170 {
171 REPORT // 426
172 Real* Mstore = store+mrc.rowcol; int i=nrows_val;
173 Real* Cstore = mrc.data;
174 // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_val; }
175 if (i) for (;;)
176 { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_val; }
177 }
178}
179
180void Matrix::RestoreCol(MatrixColX& mrc)
181{
182 REPORT
183 Real* Mstore = store+mrc.rowcol; int i=nrows_val; Real* Cstore = mrc.data;
184 // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_val; }
185 if (i) for (;;)
186 { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_val; }
187}
188
189void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); } // 1808
190
191void Matrix::NextCol(MatrixRowCol& mrc)
192{
193 REPORT // 632
194 if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
195 mrc.rowcol++;
196 if (mrc.rowcol<ncols_val)
197 {
198 if (+(mrc.cw*LoadOnEntry))
199 {
200 REPORT
201 Real* ColCopy = mrc.data;
202 Real* Mstore = store+mrc.rowcol; int i=nrows_val;
203 //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
204 if (i) for (;;)
205 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
206 }
207 }
208 else { REPORT mrc.cw -= StoreOnExit; }
209}
210
211void Matrix::NextCol(MatrixColX& mrc)
212{
213 REPORT
214 if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
215 mrc.rowcol++;
216 if (mrc.rowcol<ncols_val)
217 {
218 if (+(mrc.cw*LoadOnEntry))
219 {
220 REPORT
221 Real* ColCopy = mrc.data;
222 Real* Mstore = store+mrc.rowcol; int i=nrows_val;
223 // while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
224 if (i) for (;;)
225 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
226 }
227 }
228 else { REPORT mrc.cw -= StoreOnExit; }
229}
230
231// routines for diagonal matrix
232
233void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
234{
235 REPORT
236 mrc.skip=mrc.rowcol; mrc.storage=1;
237 mrc.data=store+mrc.skip; mrc.length=ncols_val;
238}
239
240void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
241{
242 REPORT
243 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
244 if (+(mrc.cw*StoreHere)) // should not happen
245 Throw(InternalException("DiagonalMatrix::GetCol(MatrixRowCol&)"));
246 else { REPORT mrc.data=store+mrc.skip; }
247 // not accessed
248}
249
250void DiagonalMatrix::GetCol(MatrixColX& mrc)
251{
252 REPORT
253 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
254 mrc.data = mrc.store+mrc.skip;
255 *(mrc.data)=*(store+mrc.skip);
256}
257
258void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
259 // 800
260
261void DiagonalMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
262 // not accessed
263
264void DiagonalMatrix::NextCol(MatrixColX& mrc)
265{
266 REPORT
267 if (+(mrc.cw*StoreOnExit))
268 { REPORT *(store+mrc.rowcol)=*(mrc.data); }
269 mrc.IncrDiag();
270 int t1 = +(mrc.cw*LoadOnEntry);
271 if (t1 && mrc.rowcol < ncols_val)
272 { REPORT *(mrc.data)=*(store+mrc.rowcol); }
273}
274
275// routines for upper triangular matrix
276
277void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
278{
279 REPORT
280 int row = mrc.rowcol; mrc.skip=row; mrc.length=ncols_val;
281 mrc.storage=ncols_val-row; mrc.data=store+(row*(2*ncols_val-row+1))/2;
282}
283
284
285void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
286{
287 REPORT
288 mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
289 mrc.length=nrows_val; Real* ColCopy;
290 if ( !(mrc.cw*(StoreHere+HaveStore)) )
291 {
292 REPORT // not accessed
293 ColCopy = new Real [nrows_val]; MatrixErrorNoSpace(ColCopy);
294 MONITOR_REAL_NEW("Make (UT GetCol)",nrows_val,ColCopy)
295 mrc.data = ColCopy; mrc.cw += HaveStore;
296 }
297 else { REPORT ColCopy = mrc.data; }
298 if (+(mrc.cw*LoadOnEntry))
299 {
300 REPORT
301 Real* Mstore = store+mrc.rowcol; int j = ncols_val;
302 // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
303 if (i) for (;;)
304 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += --j; }
305 }
306}
307
308void UpperTriangularMatrix::GetCol(MatrixColX& mrc)
309{
310 REPORT
311 mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
312 mrc.length=nrows_val;
313 if (+(mrc.cw*LoadOnEntry))
314 {
315 REPORT
316 Real* ColCopy = mrc.data;
317 Real* Mstore = store+mrc.rowcol; int j = ncols_val;
318 // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
319 if (i) for (;;)
320 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += --j; }
321 }
322}
323
324void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
325{
326 REPORT
327 Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols_val;
328 Real* Cstore = mrc.data;
329 // while (i--) { *Mstore = *Cstore++; Mstore += --j; }
330 if (i) for (;;)
331 { *Mstore = *Cstore++; if (!(--i)) break; Mstore += --j; }
332}
333
334void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
335 // 722
336
337// routines for lower triangular matrix
338
339void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
340{
341 REPORT
342 int row=mrc.rowcol; mrc.skip=0; mrc.storage=row+1; mrc.length=ncols_val;
343 mrc.data=store+(row*(row+1))/2;
344}
345
346void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
347{
348 REPORT
349 int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows_val;
350 int i=nrows_val-col; mrc.storage=i; Real* ColCopy;
351 if ( +(mrc.cw*(StoreHere+HaveStore)) )
352 { REPORT ColCopy = mrc.data; }
353 else
354 {
355 REPORT // not accessed
356 ColCopy = new Real [nrows_val]; MatrixErrorNoSpace(ColCopy);
357 MONITOR_REAL_NEW("Make (LT GetCol)",nrows_val,ColCopy)
358 mrc.cw += HaveStore; mrc.data = ColCopy;
359 }
360
361 if (+(mrc.cw*LoadOnEntry))
362 {
363 REPORT
364 Real* Mstore = store+(col*(col+3))/2;
365 // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
366 if (i) for (;;)
367 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
368 }
369}
370
371void LowerTriangularMatrix::GetCol(MatrixColX& mrc)
372{
373 REPORT
374 int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows_val;
375 int i=nrows_val-col; mrc.storage=i; mrc.data = mrc.store + col;
376
377 if (+(mrc.cw*LoadOnEntry))
378 {
379 REPORT Real* ColCopy = mrc.data;
380 Real* Mstore = store+(col*(col+3))/2;
381 // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
382 if (i) for (;;)
383 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
384 }
385}
386
387void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
388{
389 REPORT
390 int col=mrc.rowcol; Real* Cstore = mrc.data;
391 Real* Mstore = store+(col*(col+3))/2; int i=nrows_val-col;
392 //while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
393 if (i) for (;;)
394 { *Mstore = *Cstore++; if (!(--i)) break; Mstore += ++col; }
395}
396
397void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
398 //712
399
400// routines for symmetric matrix
401
402void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
403{
404 REPORT //571
405 mrc.skip=0; int row=mrc.rowcol; mrc.length=ncols_val;
406 if (+(mrc.cw*DirectPart))
407 { REPORT mrc.storage=row+1; mrc.data=store+(row*(row+1))/2; }
408 else
409 {
410 // do not allow StoreOnExit and !DirectPart
411 if (+(mrc.cw*StoreOnExit))
412 Throw(InternalException("SymmetricMatrix::GetRow(MatrixRowCol&)"));
413 mrc.storage=ncols_val; Real* RowCopy;
414 if (!(mrc.cw*HaveStore))
415 {
416 REPORT
417 RowCopy = new Real [ncols_val]; MatrixErrorNoSpace(RowCopy);
418 MONITOR_REAL_NEW("Make (SymGetRow)",ncols_val,RowCopy)
419 mrc.cw += HaveStore; mrc.data = RowCopy;
420 }
421 else { REPORT RowCopy = mrc.data; }
422 if (+(mrc.cw*LoadOnEntry))
423 {
424 REPORT // 544
425 Real* Mstore = store+(row*(row+1))/2; int i = row;
426 while (i--) *RowCopy++ = *Mstore++;
427 i = ncols_val-row;
428 // while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
429 if (i) for (;;)
430 { *RowCopy++ = *Mstore; if (!(--i)) break; Mstore += ++row; }
431 }
432 }
433}
434
435void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
436{
437 // do not allow StoreHere
438 if (+(mrc.cw*StoreHere))
439 Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
440
441 int col=mrc.rowcol; mrc.length=nrows_val;
442 REPORT
443 mrc.skip=0;
444 if (+(mrc.cw*DirectPart)) // actually get row ??
445 { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
446 else
447 {
448 // do not allow StoreOnExit and !DirectPart
449 if (+(mrc.cw*StoreOnExit))
450 Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
451
452 mrc.storage=ncols_val; Real* ColCopy;
453 if ( +(mrc.cw*HaveStore)) { REPORT ColCopy = mrc.data; }
454 else
455 {
456 REPORT // not accessed
457 ColCopy = new Real [ncols_val]; MatrixErrorNoSpace(ColCopy);
458 MONITOR_REAL_NEW("Make (SymGetCol)",ncols_val,ColCopy)
459 mrc.cw += HaveStore; mrc.data = ColCopy;
460 }
461 if (+(mrc.cw*LoadOnEntry))
462 {
463 REPORT
464 Real* Mstore = store+(col*(col+1))/2; int i = col;
465 while (i--) *ColCopy++ = *Mstore++;
466 i = ncols_val-col;
467 // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
468 if (i) for (;;)
469 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
470 }
471 }
472}
473
474void SymmetricMatrix::GetCol(MatrixColX& mrc)
475{
476 int col=mrc.rowcol; mrc.length=nrows_val;
477 if (+(mrc.cw*DirectPart))
478 {
479 REPORT
480 mrc.skip=col; int i=nrows_val-col; mrc.storage=i;
481 mrc.data = mrc.store+col;
482 if (+(mrc.cw*LoadOnEntry))
483 {
484 REPORT // not accessed
485 Real* ColCopy = mrc.data;
486 Real* Mstore = store+(col*(col+3))/2;
487 // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
488 if (i) for (;;)
489 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
490 }
491 }
492 else
493 {
494 REPORT
495 // do not allow StoreOnExit and !DirectPart
496 if (+(mrc.cw*StoreOnExit))
497 Throw(InternalException("SymmetricMatrix::GetCol(MatrixColX&)"));
498
499 mrc.skip=0; mrc.storage=ncols_val;
500 if (+(mrc.cw*LoadOnEntry))
501 {
502 REPORT
503 Real* ColCopy = mrc.data;
504 Real* Mstore = store+(col*(col+1))/2; int i = col;
505 while (i--) *ColCopy++ = *Mstore++;
506 i = ncols_val-col;
507 // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
508 if (i) for (;;)
509 { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
510 }
511 }
512}
513
514// Do not need RestoreRow because we do not allow !DirectPart && StoreOnExit
515
516void SymmetricMatrix::RestoreCol(MatrixColX& mrc)
517{
518 REPORT
519 // Really do restore column
520 int col=mrc.rowcol; Real* Cstore = mrc.data;
521 Real* Mstore = store+(col*(col+3))/2; int i = nrows_val-col;
522 // while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
523 if (i) for (;;)
524 { *Mstore = *Cstore++; if (!(--i)) break; Mstore+= ++col; }
525
526}
527
528// routines for row vector
529
530void RowVector::GetCol(MatrixRowCol& mrc)
531{
532 REPORT
533 // do not allow StoreHere
534 if (+(mrc.cw*StoreHere))
535 Throw(InternalException("RowVector::GetCol(MatrixRowCol&)"));
536
537 mrc.skip=0; mrc.storage=1; mrc.length=nrows_val;
538 mrc.data = store+mrc.rowcol;
539
540}
541
542void RowVector::GetCol(MatrixColX& mrc)
543{
544 REPORT
545 mrc.skip=0; mrc.storage=1; mrc.length=nrows_val;
546 if (mrc.cw >= LoadOnEntry)
547 { REPORT *(mrc.data) = *(store+mrc.rowcol); }
548
549}
550
551void RowVector::NextCol(MatrixRowCol& mrc)
552{ REPORT mrc.rowcol++; mrc.data++; }
553
554void RowVector::NextCol(MatrixColX& mrc)
555{
556 if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.data); }
557
558 mrc.rowcol++;
559 if (mrc.rowcol < ncols_val)
560 {
561 if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.data)=*(store+mrc.rowcol); }
562 }
563 else { REPORT mrc.cw -= StoreOnExit; }
564}
565
566void RowVector::RestoreCol(MatrixColX& mrc)
567 { REPORT *(store+mrc.rowcol)=*(mrc.data); } // not accessed
568
569
570// routines for band matrices
571
572void BandMatrix::GetRow(MatrixRowCol& mrc)
573{
574 REPORT
575 int r = mrc.rowcol; int w = lower_val+1+upper_val; mrc.length=ncols_val;
576 int s = r-lower_val;
577 if (s<0) { mrc.data = store+(r*w-s); w += s; s = 0; }
578 else mrc.data = store+r*w;
579 mrc.skip = s; s += w-ncols_val; if (s>0) w -= s; mrc.storage = w;
580}
581
582// should make special versions of this for upper and lower band matrices
583
584void BandMatrix::NextRow(MatrixRowCol& mrc)
585{
586 REPORT
587 int r = ++mrc.rowcol;
588 if (r<=lower_val) { mrc.storage++; mrc.data += lower_val+upper_val; }
589 else { mrc.skip++; mrc.data += lower_val+upper_val+1; }
590 if (r>=ncols_val-upper_val) mrc.storage--;
591}
592
593void BandMatrix::GetCol(MatrixRowCol& mrc)
594{
595 REPORT
596 int c = mrc.rowcol; int n = lower_val+upper_val; int w = n+1;
597 mrc.length=nrows_val; Real* ColCopy;
598 int b; int s = c-upper_val;
599 if (s<=0) { w += s; s = 0; b = c+lower_val; } else b = s*w+n;
600 mrc.skip = s; s += w-nrows_val; if (s>0) w -= s; mrc.storage = w;
601 if ( +(mrc.cw*(StoreHere+HaveStore)) )
602 { REPORT ColCopy = mrc.data; }
603 else
604 {
605 REPORT
606 ColCopy = new Real [n+1]; MatrixErrorNoSpace(ColCopy);
607 MONITOR_REAL_NEW("Make (BMGetCol)",n+1,ColCopy)
608 mrc.cw += HaveStore; mrc.data = ColCopy;
609 }
610
611 if (+(mrc.cw*LoadOnEntry))
612 {
613 REPORT
614 Real* Mstore = store+b;
615 // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
616 if (w) for (;;)
617 { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
618 }
619}
620
621void BandMatrix::GetCol(MatrixColX& mrc)
622{
623 REPORT
624 int c = mrc.rowcol; int n = lower_val+upper_val; int w = n+1;
625 mrc.length=nrows_val; int b; int s = c-upper_val;
626 if (s<=0) { w += s; s = 0; b = c+lower_val; } else b = s*w+n;
627 mrc.skip = s; s += w-nrows_val; if (s>0) w -= s; mrc.storage = w;
628 mrc.data = mrc.store+mrc.skip;
629
630 if (+(mrc.cw*LoadOnEntry))
631 {
632 REPORT
633 Real* ColCopy = mrc.data; Real* Mstore = store+b;
634 // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
635 if (w) for (;;)
636 { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
637 }
638}
639
640void BandMatrix::RestoreCol(MatrixRowCol& mrc)
641{
642 REPORT
643 int c = mrc.rowcol; int n = lower_val+upper_val; int s = c-upper_val;
644 Real* Mstore = store + ((s<=0) ? c+lower_val : s*n+s+n);
645 Real* Cstore = mrc.data;
646 int w = mrc.storage;
647 // while (w--) { *Mstore = *Cstore++; Mstore += n; }
648 if (w) for (;;)
649 { *Mstore = *Cstore++; if (!(--w)) break; Mstore += n; }
650}
651
652// routines for symmetric band matrix
653
654void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
655{
656 REPORT
657 int r=mrc.rowcol; int s = r-lower_val; int w1 = lower_val+1; int o = r*w1;
658 mrc.length = ncols_val;
659 if (s<0) { w1 += s; o -= s; s = 0; }
660 mrc.skip = s;
661
662 if (+(mrc.cw*DirectPart))
663 { REPORT mrc.data = store+o; mrc.storage = w1; }
664 else
665 {
666 // do not allow StoreOnExit and !DirectPart
667 if (+(mrc.cw*StoreOnExit))
668 Throw(InternalException("SymmetricBandMatrix::GetRow(MatrixRowCol&)"));
669 int w = w1+lower_val; s += w-ncols_val; Real* RowCopy;
670 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
671 if (!(mrc.cw*HaveStore))
672 {
673 REPORT
674 RowCopy = new Real [2*lower_val+1]; MatrixErrorNoSpace(RowCopy);
675 MONITOR_REAL_NEW("Make (SmBGetRow)",2*lower_val+1,RowCopy)
676 mrc.cw += HaveStore; mrc.data = RowCopy;
677 }
678 else { REPORT RowCopy = mrc.data; }
679
680 if (+(mrc.cw*LoadOnEntry) && ncols_val > 0)
681 {
682 REPORT
683 Real* Mstore = store+o;
684 while (w1--) *RowCopy++ = *Mstore++;
685 Mstore--;
686 while (w2--) { Mstore += lower_val; *RowCopy++ = *Mstore; }
687 }
688 }
689}
690
691void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
692{
693 // do not allow StoreHere
694 if (+(mrc.cw*StoreHere))
695 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
696
697 int c=mrc.rowcol; int w1 = lower_val+1; mrc.length=nrows_val;
698 REPORT
699 int s = c-lower_val; int o = c*w1;
700 if (s<0) { w1 += s; o -= s; s = 0; }
701 mrc.skip = s;
702
703 if (+(mrc.cw*DirectPart))
704 { REPORT mrc.data = store+o; mrc.storage = w1; }
705 else
706 {
707 // do not allow StoreOnExit and !DirectPart
708 if (+(mrc.cw*StoreOnExit))
709 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
710 int w = w1+lower_val; s += w-ncols_val; Real* ColCopy;
711 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
712
713 if ( +(mrc.cw*HaveStore) ) { REPORT ColCopy = mrc.data; }
714 else
715 {
716 REPORT ColCopy = new Real [2*lower_val+1]; MatrixErrorNoSpace(ColCopy);
717 MONITOR_REAL_NEW("Make (SmBGetCol)",2*lower_val+1,ColCopy)
718 mrc.cw += HaveStore; mrc.data = ColCopy;
719 }
720
721 if (+(mrc.cw*LoadOnEntry))
722 {
723 REPORT
724 Real* Mstore = store+o;
725 while (w1--) *ColCopy++ = *Mstore++;
726 Mstore--;
727 while (w2--) { Mstore += lower_val; *ColCopy++ = *Mstore; }
728 }
729 }
730}
731
732void SymmetricBandMatrix::GetCol(MatrixColX& mrc)
733{
734 int c=mrc.rowcol; int w1 = lower_val+1; mrc.length=nrows_val;
735 if (+(mrc.cw*DirectPart))
736 {
737 REPORT
738 int b = c*w1+lower_val;
739 mrc.skip = c; c += w1-nrows_val; w1 -= c; mrc.storage = w1;
740 Real* ColCopy = mrc.data = mrc.store+mrc.skip;
741 if (+(mrc.cw*LoadOnEntry))
742 {
743 REPORT
744 Real* Mstore = store+b;
745 // while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
746 if (w1) for (;;)
747 { *ColCopy++ = *Mstore; if (!(--w1)) break; Mstore += lower_val; }
748 }
749 }
750 else
751 {
752 REPORT
753 // do not allow StoreOnExit and !DirectPart
754 if (+(mrc.cw*StoreOnExit))
755 Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixColX&)"));
756 int s = c-lower_val; int o = c*w1;
757 if (s<0) { w1 += s; o -= s; s = 0; }
758 mrc.skip = s;
759
760 int w = w1+lower_val; s += w-ncols_val;
761 if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
762
763 Real* ColCopy = mrc.data = mrc.store+mrc.skip;
764
765 if (+(mrc.cw*LoadOnEntry))
766 {
767 REPORT
768 Real* Mstore = store+o;
769 while (w1--) *ColCopy++ = *Mstore++;
770 Mstore--;
771 while (w2--) { Mstore += lower_val; *ColCopy++ = *Mstore; }
772 }
773
774 }
775}
776
777void SymmetricBandMatrix::RestoreCol(MatrixColX& mrc)
778{
779 REPORT
780 int c = mrc.rowcol;
781 Real* Mstore = store + c*lower_val+c+lower_val;
782 Real* Cstore = mrc.data; int w = mrc.storage;
783 // while (w--) { *Mstore = *Cstore++; Mstore += lower_val; }
784 if (w) for (;;)
785 { *Mstore = *Cstore++; if (!(--w)) break; Mstore += lower_val; }
786}
787
788// routines for identity matrix
789
790void IdentityMatrix::GetRow(MatrixRowCol& mrc)
791{
792 REPORT
793 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.data=store; mrc.length=ncols_val;
794}
795
796void IdentityMatrix::GetCol(MatrixRowCol& mrc)
797{
798 REPORT
799 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
800 if (+(mrc.cw*StoreHere)) // should not happen
801 Throw(InternalException("IdentityMatrix::GetCol(MatrixRowCol&)"));
802 else { REPORT mrc.data=store; }
803}
804
805void IdentityMatrix::GetCol(MatrixColX& mrc)
806{
807 REPORT
808 mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
809 mrc.data = mrc.store+mrc.skip; *(mrc.data)=*store;
810}
811
812void IdentityMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrId(); }
813
814void IdentityMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrId(); }
815
816void IdentityMatrix::NextCol(MatrixColX& mrc)
817{
818 REPORT
819 if (+(mrc.cw*StoreOnExit)) { REPORT *store=*(mrc.data); }
820 mrc.IncrDiag(); // must increase mrc.data so need IncrDiag
821 int t1 = +(mrc.cw*LoadOnEntry);
822 if (t1 && mrc.rowcol < ncols_val) { REPORT *(mrc.data)=*store; }
823}
824
825
826
827
828// *************************** destructors *******************************
829
830MatrixRowCol::~MatrixRowCol()
831{
832 if (+(cw*HaveStore))
833 {
834 MONITOR_REAL_DELETE("Free (RowCol)",-1,data) // do not know length
835 delete [] data;
836 }
837}
838
839MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }
840
841MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
842
843MatrixColX::~MatrixColX() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
844
845#ifdef use_namespace
846}
847#endif
848
849///@}
Note: See TracBrowser for help on using the repository browser.