Context Navigation

newmat3.cpp@ 9448

Last change on this file since 9448 was 9434, checked in by stuerze, 4 years ago
newmat is reverted to its stable version 10 because version 11beta was not working with mac compilers
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
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(2ncols_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.cwStoreOnExit)) { REPORT (store+mrc.rowcol)=*(mrc.data); }
557
558	mrc.rowcol++;
559	if (mrc.rowcol < ncols_val)
560	{
561	if (+(mrc.cwLoadOnEntry)) { 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.cwStoreOnExit)) { 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 (+(cwStoreOnExit)) gm->RestoreRow(this); }
840
841	MatrixCol::~MatrixCol() { if (+(cwStoreOnExit)) gm->RestoreCol(this); }
842
843	MatrixColX::~MatrixColX() { if (+(cwStoreOnExit)) gm->RestoreCol(this); }
844
845	#ifdef use_namespace
846	}
847	#endif
848
849	///@}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: ntrip/trunk/BNC/newmat/newmat3.cpp@ 9448

Download in other formats: