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 |
---|
17 | namespace 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 | |
---|
85 | void 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 | |
---|
94 | void 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 | |
---|
103 | void 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 | |
---|
115 | void 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 | |
---|
123 | void 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 | |
---|
151 | void 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 | |
---|
165 | void 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 | |
---|
180 | void 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 | |
---|
189 | void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); } // 1808 |
---|
190 | |
---|
191 | void 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 | |
---|
211 | void 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 | |
---|
233 | void 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 | |
---|
240 | void 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 | |
---|
250 | void 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 | |
---|
258 | void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); } |
---|
259 | // 800 |
---|
260 | |
---|
261 | void DiagonalMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); } |
---|
262 | // not accessed |
---|
263 | |
---|
264 | void 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 | |
---|
277 | void 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 | |
---|
285 | void 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 | |
---|
308 | void 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 | |
---|
324 | void 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 | |
---|
334 | void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); } |
---|
335 | // 722 |
---|
336 | |
---|
337 | // routines for lower triangular matrix |
---|
338 | |
---|
339 | void 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 | |
---|
346 | void 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 | |
---|
371 | void 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 | |
---|
387 | void 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 | |
---|
397 | void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); } |
---|
398 | //712 |
---|
399 | |
---|
400 | // routines for symmetric matrix |
---|
401 | |
---|
402 | void 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 | |
---|
435 | void 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 | |
---|
474 | void 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 | |
---|
516 | void 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 | |
---|
530 | void 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 | |
---|
542 | void 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 | |
---|
551 | void RowVector::NextCol(MatrixRowCol& mrc) |
---|
552 | { REPORT mrc.rowcol++; mrc.data++; } |
---|
553 | |
---|
554 | void 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 | |
---|
566 | void RowVector::RestoreCol(MatrixColX& mrc) |
---|
567 | { REPORT *(store+mrc.rowcol)=*(mrc.data); } // not accessed |
---|
568 | |
---|
569 | |
---|
570 | // routines for band matrices |
---|
571 | |
---|
572 | void 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 | |
---|
584 | void 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 | |
---|
593 | void 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 | |
---|
621 | void 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 | |
---|
640 | void 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 | |
---|
654 | void 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 | |
---|
691 | void 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 | |
---|
732 | void 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 | |
---|
777 | void 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 | |
---|
790 | void IdentityMatrix::GetRow(MatrixRowCol& mrc) |
---|
791 | { |
---|
792 | REPORT |
---|
793 | mrc.skip=mrc.rowcol; mrc.storage=1; mrc.data=store; mrc.length=ncols_val; |
---|
794 | } |
---|
795 | |
---|
796 | void 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 | |
---|
805 | void 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 | |
---|
812 | void IdentityMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrId(); } |
---|
813 | |
---|
814 | void IdentityMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrId(); } |
---|
815 | |
---|
816 | void 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 | |
---|
830 | MatrixRowCol::~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 | |
---|
839 | MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); } |
---|
840 | |
---|
841 | MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); } |
---|
842 | |
---|
843 | MatrixColX::~MatrixColX() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); } |
---|
844 | |
---|
845 | #ifdef use_namespace |
---|
846 | } |
---|
847 | #endif |
---|
848 | |
---|
849 | ///@} |
---|