1 /// \ingroup newmat
2 ///@{
3
4 /// \file newmat2.cpp
5 /// Matrix row and column operations.
6 /// The operations on individual rows and columns used to carry out matrix
7 /// add, multiply etc.
8
9
10 // Copyright (C) 1991,2,3,4: R B Davies
11
12 #define WANT_MATH
13
14 #include "include.h"
15
16 #include "newmat.h"
17 #include "newmatrc.h"
18
19 #ifdef use_namespace
20 namespace NEWMAT {
21 #endif
22
23
24 #ifdef DO_REPORT
25 #define REPORT { static ExeCounter ExeCount(__LINE__,2); ++ExeCount; }
26 #else
27 #define REPORT {}
28 #endif
29
30 //#define MONITOR(what,storage,store) { cout << what << " " << storage << " at " << (long)store << "\n"; }
31
32 #define MONITOR(what,store,storage) {}
33
34 /************************** Matrix Row/Col functions ************************/
35
Add(const MatrixRowCol & mrc)36 void MatrixRowCol::Add(const MatrixRowCol& mrc)
37 {
38 // THIS += mrc
39 REPORT
40 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
41 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
42 if (l<=0) return;
43 Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
44 while (l--) *elx++ += *el++;
45 }
46
AddScaled(const MatrixRowCol & mrc,Real x)47 void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x)
48 {
49 REPORT
50 // THIS += (mrc * x)
51 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
52 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
53 if (l<=0) return;
54 Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
55 while (l--) *elx++ += *el++ * x;
56 }
57
Sub(const MatrixRowCol & mrc)58 void MatrixRowCol::Sub(const MatrixRowCol& mrc)
59 {
60 REPORT
61 // THIS -= mrc
62 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
63 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
64 if (l<=0) return;
65 Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
66 while (l--) *elx++ -= *el++;
67 }
68
Inject(const MatrixRowCol & mrc)69 void MatrixRowCol::Inject(const MatrixRowCol& mrc)
70 // copy stored elements only
71 {
72 REPORT
73 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
74 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
75 if (l<=0) return;
76 Real* elx=data+(f-skip); Real* ely=mrc.data+(f-mrc.skip);
77 while (l--) *elx++ = *ely++;
78 }
79
DotProd(const MatrixRowCol & mrc1,const MatrixRowCol & mrc2)80 Real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
81 {
82 REPORT // not accessed
83 int f = mrc1.skip; int f2 = mrc2.skip;
84 int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
85 if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
86 if (l<=0) return 0.0;
87
88 Real* el1=mrc1.data+(f-mrc1.skip); Real* el2=mrc2.data+(f-mrc2.skip);
89 Real sum = 0.0;
90 while (l--) sum += *el1++ * *el2++;
91 return sum;
92 }
93
Add(const MatrixRowCol & mrc1,const MatrixRowCol & mrc2)94 void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
95 {
96 // THIS = mrc1 + mrc2
97 int f = skip; int l = skip + storage;
98 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
99 if (f1<f) f1=f; if (l1>l) l1=l;
100 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
101 if (f2<f) f2=f; if (l2>l) l2=l;
102 Real* el = data + (f-skip);
103 Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
104 if (f1<f2)
105 {
106 int i = f1-f; while (i--) *el++ = 0.0;
107 if (l1<=f2) // disjoint
108 {
109 REPORT // not accessed
110 i = l1-f1; while (i--) *el++ = *el1++;
111 i = f2-l1; while (i--) *el++ = 0.0;
112 i = l2-f2; while (i--) *el++ = *el2++;
113 i = l-l2; while (i--) *el++ = 0.0;
114 }
115 else
116 {
117 i = f2-f1; while (i--) *el++ = *el1++;
118 if (l1<=l2)
119 {
120 REPORT
121 i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
122 i = l2-l1; while (i--) *el++ = *el2++;
123 i = l-l2; while (i--) *el++ = 0.0;
124 }
125 else
126 {
127 REPORT
128 i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
129 i = l1-l2; while (i--) *el++ = *el1++;
130 i = l-l1; while (i--) *el++ = 0.0;
131 }
132 }
133 }
134 else
135 {
136 int i = f2-f; while (i--) *el++ = 0.0;
137 if (l2<=f1) // disjoint
138 {
139 REPORT // not accessed
140 i = l2-f2; while (i--) *el++ = *el2++;
141 i = f1-l2; while (i--) *el++ = 0.0;
142 i = l1-f1; while (i--) *el++ = *el1++;
143 i = l-l1; while (i--) *el++ = 0.0;
144 }
145 else
146 {
147 i = f1-f2; while (i--) *el++ = *el2++;
148 if (l2<=l1)
149 {
150 REPORT
151 i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
152 i = l1-l2; while (i--) *el++ = *el1++;
153 i = l-l1; while (i--) *el++ = 0.0;
154 }
155 else
156 {
157 REPORT
158 i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
159 i = l2-l1; while (i--) *el++ = *el2++;
160 i = l-l2; while (i--) *el++ = 0.0;
161 }
162 }
163 }
164 }
165
Sub(const MatrixRowCol & mrc1,const MatrixRowCol & mrc2)166 void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
167 {
168 // THIS = mrc1 - mrc2
169 int f = skip; int l = skip + storage;
170 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
171 if (f1<f) f1=f; if (l1>l) l1=l;
172 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
173 if (f2<f) f2=f; if (l2>l) l2=l;
174 Real* el = data + (f-skip);
175 Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
176 if (f1<f2)
177 {
178 int i = f1-f; while (i--) *el++ = 0.0;
179 if (l1<=f2) // disjoint
180 {
181 REPORT // not accessed
182 i = l1-f1; while (i--) *el++ = *el1++;
183 i = f2-l1; while (i--) *el++ = 0.0;
184 i = l2-f2; while (i--) *el++ = - *el2++;
185 i = l-l2; while (i--) *el++ = 0.0;
186 }
187 else
188 {
189 i = f2-f1; while (i--) *el++ = *el1++;
190 if (l1<=l2)
191 {
192 REPORT
193 i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
194 i = l2-l1; while (i--) *el++ = - *el2++;
195 i = l-l2; while (i--) *el++ = 0.0;
196 }
197 else
198 {
199 REPORT
200 i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
201 i = l1-l2; while (i--) *el++ = *el1++;
202 i = l-l1; while (i--) *el++ = 0.0;
203 }
204 }
205 }
206 else
207 {
208 int i = f2-f; while (i--) *el++ = 0.0;
209 if (l2<=f1) // disjoint
210 {
211 REPORT // not accessed
212 i = l2-f2; while (i--) *el++ = - *el2++;
213 i = f1-l2; while (i--) *el++ = 0.0;
214 i = l1-f1; while (i--) *el++ = *el1++;
215 i = l-l1; while (i--) *el++ = 0.0;
216 }
217 else
218 {
219 i = f1-f2; while (i--) *el++ = - *el2++;
220 if (l2<=l1)
221 {
222 REPORT
223 i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
224 i = l1-l2; while (i--) *el++ = *el1++;
225 i = l-l1; while (i--) *el++ = 0.0;
226 }
227 else
228 {
229 REPORT
230 i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
231 i = l2-l1; while (i--) *el++ = - *el2++;
232 i = l-l2; while (i--) *el++ = 0.0;
233 }
234 }
235 }
236 }
237
238
Add(const MatrixRowCol & mrc1,Real x)239 void MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x)
240 {
241 // THIS = mrc1 + x
242 REPORT
243 if (!storage) return;
244 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
245 if (f < skip) { f = skip; if (l < f) l = f; }
246 if (l > lx) { l = lx; if (f > lx) f = lx; }
247
248 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
249
250 int l1 = f-skip; while (l1--) *elx++ = x;
251 l1 = l-f; while (l1--) *elx++ = *ely++ + x;
252 lx -= l; while (lx--) *elx++ = x;
253 }
254
NegAdd(const MatrixRowCol & mrc1,Real x)255 void MatrixRowCol::NegAdd(const MatrixRowCol& mrc1, Real x)
256 {
257 // THIS = x - mrc1
258 REPORT
259 if (!storage) return;
260 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
261 if (f < skip) { f = skip; if (l < f) l = f; }
262 if (l > lx) { l = lx; if (f > lx) f = lx; }
263
264 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
265
266 int l1 = f-skip; while (l1--) *elx++ = x;
267 l1 = l-f; while (l1--) *elx++ = x - *ely++;
268 lx -= l; while (lx--) *elx++ = x;
269 }
270
RevSub(const MatrixRowCol & mrc1)271 void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
272 {
273 // THIS = mrc - THIS
274 REPORT
275 if (!storage) return;
276 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
277 if (f < skip) { f = skip; if (l < f) l = f; }
278 if (l > lx) { l = lx; if (f > lx) f = lx; }
279
280 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
281
282 int l1 = f-skip; while (l1--) { *elx = - *elx; elx++; }
283 l1 = l-f; while (l1--) { *elx = *ely++ - *elx; elx++; }
284 lx -= l; while (lx--) { *elx = - *elx; elx++; }
285 }
286
ConCat(const MatrixRowCol & mrc1,const MatrixRowCol & mrc2)287 void MatrixRowCol::ConCat(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
288 {
289 // THIS = mrc1 | mrc2
290 REPORT
291 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage; int lx = skip + storage;
292 if (f1 < skip) { f1 = skip; if (l1 < f1) l1 = f1; }
293 if (l1 > lx) { l1 = lx; if (f1 > lx) f1 = lx; }
294
295 Real* elx = data;
296
297 int i = f1-skip; while (i--) *elx++ =0.0;
298 i = l1-f1;
299 if (i) // in case f1 would take ely out of range
300 { Real* ely = mrc1.data+(f1-mrc1.skip); while (i--) *elx++ = *ely++; }
301
302 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage; i = mrc1.length;
303 int skipx = l1 - i; lx -= i; // addresses rel to second seg, maybe -ve
304 if (f2 < skipx) { f2 = skipx; if (l2 < f2) l2 = f2; }
305 if (l2 > lx) { l2 = lx; if (f2 > lx) f2 = lx; }
306
307 i = f2-skipx; while (i--) *elx++ = 0.0;
308 i = l2-f2;
309 if (i) // in case f2 would take ely out of range
310 { Real* ely = mrc2.data+(f2-mrc2.skip); while (i--) *elx++ = *ely++; }
311 lx -= l2; while (lx--) *elx++ = 0.0;
312 }
313
Multiply(const MatrixRowCol & mrc1)314 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1)
315 // element by element multiply into
316 {
317 REPORT
318 if (!storage) return;
319 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
320 if (f < skip) { f = skip; if (l < f) l = f; }
321 if (l > lx) { l = lx; if (f > lx) f = lx; }
322
323 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
324
325 int l1 = f-skip; while (l1--) *elx++ = 0;
326 l1 = l-f; while (l1--) *elx++ *= *ely++;
327 lx -= l; while (lx--) *elx++ = 0;
328 }
329
Multiply(const MatrixRowCol & mrc1,const MatrixRowCol & mrc2)330 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
331 // element by element multiply
332 {
333 int f = skip; int l = skip + storage;
334 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
335 if (f1<f) f1=f; if (l1>l) l1=l;
336 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
337 if (f2<f) f2=f; if (l2>l) l2=l;
338 Real* el = data + (f-skip); int i;
339 if (f1<f2) f1 = f2; if (l1>l2) l1 = l2;
340 if (l1<=f1) { REPORT i = l-f; while (i--) *el++ = 0.0; } // disjoint
341 else
342 {
343 REPORT
344 Real* el1 = mrc1.data+(f1-mrc1.skip);
345 Real* el2 = mrc2.data+(f1-mrc2.skip);
346 i = f1-f ; while (i--) *el++ = 0.0;
347 i = l1-f1; while (i--) *el++ = *el1++ * *el2++;
348 i = l-l1; while (i--) *el++ = 0.0;
349 }
350 }
351
KP(const MatrixRowCol & mrc1,const MatrixRowCol & mrc2)352 void MatrixRowCol::KP(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
353 // row for Kronecker product
354 {
355 int f = skip; int s = storage; Real* el = data; int i;
356
357 i = mrc1.skip * mrc2.length;
358 if (i > f)
359 {
360 i -= f; f = 0; if (i > s) { i = s; s = 0; } else s -= i;
361 while (i--) *el++ = 0.0;
362 if (s == 0) return;
363 }
364 else f -= i;
365
366 i = mrc1.storage; Real* el1 = mrc1.data;
367 int mrc2_skip = mrc2.skip; int mrc2_storage = mrc2.storage;
368 int mrc2_length = mrc2.length;
369 int mrc2_remain = mrc2_length - mrc2_skip - mrc2_storage;
370 while (i--)
371 {
372 int j; Real* el2 = mrc2.data; Real vel1 = *el1;
373 if (f == 0 && mrc2_length <= s)
374 {
375 j = mrc2_skip; s -= j; while (j--) *el++ = 0.0;
376 j = mrc2_storage; s -= j; while (j--) *el++ = vel1 * *el2++;
377 j = mrc2_remain; s -= j; while (j--) *el++ = 0.0;
378 }
379 else if (f >= mrc2_length) f -= mrc2_length;
380 else
381 {
382 j = mrc2_skip;
383 if (j > f)
384 {
385 j -= f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
386 while (j--) *el++ = 0.0;
387 }
388 else f -= j;
389
390 j = mrc2_storage;
391 if (j > f)
392 {
393 j -= f; el2 += f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
394 while (j--) *el++ = vel1 * *el2++;
395 }
396 else f -= j;
397
398 j = mrc2_remain;
399 if (j > f)
400 {
401 j -= f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
402 while (j--) *el++ = 0.0;
403 }
404 else f -= j;
405 }
406 if (s == 0) return;
407 ++el1;
408 }
409
410 i = (mrc1.length - mrc1.skip - mrc1.storage) * mrc2.length;
411 if (i > f)
412 {
413 i -= f; if (i > s) i = s;
414 while (i--) *el++ = 0.0;
415 }
416 }
417
418
Copy(const MatrixRowCol & mrc1)419 void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
420 {
421 // THIS = mrc1
422 REPORT
423 if (!storage) return;
424 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
425 if (f < skip) { f = skip; if (l < f) l = f; }
426 if (l > lx) { l = lx; if (f > lx) f = lx; }
427
428 Real* elx = data; Real* ely = 0;
429
430 if (l-f) ely = mrc1.data+(f-mrc1.skip);
431
432 int l1 = f-skip; while (l1--) *elx++ = 0.0;
433 l1 = l-f; while (l1--) *elx++ = *ely++;
434 lx -= l; while (lx--) *elx++ = 0.0;
435 }
436
CopyCheck(const MatrixRowCol & mrc1)437 void MatrixRowCol::CopyCheck(const MatrixRowCol& mrc1)
438 // Throw an exception if this would lead to a loss of data
439 {
440 REPORT
441 if (!storage) return;
442 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
443 if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
444
445 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
446
447 int l1 = f-skip; while (l1--) *elx++ = 0.0;
448 l1 = l-f; while (l1--) *elx++ = *ely++;
449 lx -= l; while (lx--) *elx++ = 0.0;
450 }
451
Check(const MatrixRowCol & mrc1)452 void MatrixRowCol::Check(const MatrixRowCol& mrc1)
453 // Throw an exception if +=, -=, copy etc would lead to a loss of data
454 {
455 REPORT
456 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
457 if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
458 }
459
Check()460 void MatrixRowCol::Check()
461 // Throw an exception if +=, -= of constant would lead to a loss of data
462 // that is: check full row is present
463 // may not be appropriate for symmetric matrices
464 {
465 REPORT
466 if (skip!=0 || storage!=length)
467 Throw(ProgramException("Illegal Conversion"));
468 }
469
Negate(const MatrixRowCol & mrc1)470 void MatrixRowCol::Negate(const MatrixRowCol& mrc1)
471 {
472 // THIS = -mrc1
473 REPORT
474 if (!storage) return;
475 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
476 if (f < skip) { f = skip; if (l < f) l = f; }
477 if (l > lx) { l = lx; if (f > lx) f = lx; }
478
479 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
480
481 int l1 = f-skip; while (l1--) *elx++ = 0.0;
482 l1 = l-f; while (l1--) *elx++ = - *ely++;
483 lx -= l; while (lx--) *elx++ = 0.0;
484 }
485
Multiply(const MatrixRowCol & mrc1,Real s)486 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s)
487 {
488 // THIS = mrc1 * s
489 REPORT
490 if (!storage) return;
491 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
492 if (f < skip) { f = skip; if (l < f) l = f; }
493 if (l > lx) { l = lx; if (f > lx) f = lx; }
494
495 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
496
497 int l1 = f-skip; while (l1--) *elx++ = 0.0;
498 l1 = l-f; while (l1--) *elx++ = *ely++ * s;
499 lx -= l; while (lx--) *elx++ = 0.0;
500 }
501
Solver(MatrixColX & mrc,const MatrixColX & mrc1)502 void DiagonalMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1)
503 {
504 // mrc = mrc / mrc1 (elementwise)
505 REPORT
506 int f = mrc1.skip; int f0 = mrc.skip;
507 int l = f + mrc1.storage; int lx = f0 + mrc.storage;
508 if (f < f0) { f = f0; if (l < f) l = f; }
509 if (l > lx) { l = lx; if (f > lx) f = lx; }
510
511 Real* elx = mrc.data; Real* eld = store+f;
512
513 int l1 = f-f0; while (l1--) *elx++ = 0.0;
514 l1 = l-f; while (l1--) *elx++ /= *eld++;
515 lx -= l; while (lx--) *elx++ = 0.0;
516 // Solver makes sure input and output point to same memory
517 }
518
Solver(MatrixColX & mrc,const MatrixColX & mrc1)519 void IdentityMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1)
520 {
521 // mrc = mrc / mrc1 (elementwise)
522 REPORT
523 int f = mrc1.skip; int f0 = mrc.skip;
524 int l = f + mrc1.storage; int lx = f0 + mrc.storage;
525 if (f < f0) { f = f0; if (l < f) l = f; }
526 if (l > lx) { l = lx; if (f > lx) f = lx; }
527
528 Real* elx = mrc.data; Real eldv = *store;
529
530 int l1 = f-f0; while (l1--) *elx++ = 0.0;
531 l1 = l-f; while (l1--) *elx++ /= eldv;
532 lx -= l; while (lx--) *elx++ = 0.0;
533 // Solver makes sure input and output point to same memory
534 }
535
Copy(const double * & r)536 void MatrixRowCol::Copy(const double*& r)
537 {
538 // THIS = *r
539 REPORT
540 Real* elx = data; const double* ely = r+skip; r += length;
541 int l = storage; while (l--) *elx++ = (Real)*ely++;
542 }
543
Copy(const float * & r)544 void MatrixRowCol::Copy(const float*& r)
545 {
546 // THIS = *r
547 REPORT
548 Real* elx = data; const float* ely = r+skip; r += length;
549 int l = storage; while (l--) *elx++ = (Real)*ely++;
550 }
551
Copy(const int * & r)552 void MatrixRowCol::Copy(const int*& r)
553 {
554 // THIS = *r
555 REPORT
556 Real* elx = data; const int* ely = r+skip; r += length;
557 int l = storage; while (l--) *elx++ = (Real)*ely++;
558 }
559
Copy(Real r)560 void MatrixRowCol::Copy(Real r)
561 {
562 // THIS = r
563 REPORT Real* elx = data; int l = storage; while (l--) *elx++ = r;
564 }
565
Zero()566 void MatrixRowCol::Zero()
567 {
568 // THIS = 0
569 REPORT Real* elx = data; int l = storage; while (l--) *elx++ = 0;
570 }
571
Multiply(Real r)572 void MatrixRowCol::Multiply(Real r)
573 {
574 // THIS *= r
575 REPORT Real* elx = data; int l = storage; while (l--) *elx++ *= r;
576 }
577
Add(Real r)578 void MatrixRowCol::Add(Real r)
579 {
580 // THIS += r
581 REPORT
582 Real* elx = data; int l = storage; while (l--) *elx++ += r;
583 }
584
SumAbsoluteValue()585 Real MatrixRowCol::SumAbsoluteValue()
586 {
587 REPORT
588 Real sum = 0.0; Real* elx = data; int l = storage;
589 while (l--) sum += fabs(*elx++);
590 return sum;
591 }
592
593 // max absolute value of r and elements of row/col
594 // we use <= or >= in all of these so we are sure of getting
595 // r reset at least once.
MaximumAbsoluteValue1(Real r,int & i)596 Real MatrixRowCol::MaximumAbsoluteValue1(Real r, int& i)
597 {
598 REPORT
599 Real* elx = data; int l = storage; int li = -1;
600 while (l--) { Real f = fabs(*elx++); if (r <= f) { r = f; li = l; } }
601 i = (li >= 0) ? storage - li + skip : 0;
602 return r;
603 }
604
605 // min absolute value of r and elements of row/col
MinimumAbsoluteValue1(Real r,int & i)606 Real MatrixRowCol::MinimumAbsoluteValue1(Real r, int& i)
607 {
608 REPORT
609 Real* elx = data; int l = storage; int li = -1;
610 while (l--) { Real f = fabs(*elx++); if (r >= f) { r = f; li = l; } }
611 i = (li >= 0) ? storage - li + skip : 0;
612 return r;
613 }
614
615 // max value of r and elements of row/col
Maximum1(Real r,int & i)616 Real MatrixRowCol::Maximum1(Real r, int& i)
617 {
618 REPORT
619 Real* elx = data; int l = storage; int li = -1;
620 while (l--) { Real f = *elx++; if (r <= f) { r = f; li = l; } }
621 i = (li >= 0) ? storage - li + skip : 0;
622 return r;
623 }
624
625 // min value of r and elements of row/col
Minimum1(Real r,int & i)626 Real MatrixRowCol::Minimum1(Real r, int& i)
627 {
628 REPORT
629 Real* elx = data; int l = storage; int li = -1;
630 while (l--) { Real f = *elx++; if (r >= f) { r = f; li = l; } }
631 i = (li >= 0) ? storage - li + skip : 0;
632 return r;
633 }
634
Sum()635 Real MatrixRowCol::Sum()
636 {
637 REPORT
638 Real sum = 0.0; Real* elx = data; int l = storage;
639 while (l--) sum += *elx++;
640 return sum;
641 }
642
SubRowCol(MatrixRowCol & mrc,int skip1,int l1) const643 void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
644 {
645 mrc.length = l1; int d = skip - skip1;
646 if (d<0) { mrc.skip = 0; mrc.data = data - d; }
647 else { mrc.skip = d; mrc.data = data; }
648 d = skip + storage - skip1;
649 d = ((l1 < d) ? l1 : d) - mrc.skip; mrc.storage = (d < 0) ? 0 : d;
650 mrc.cw = 0;
651 }
652
653 #ifdef use_namespace
654 }
655 #endif
656
657
658 ///@}
659