source: ntrip/trunk/BNC/newmat/newmat2.cpp@ 4636

Last change on this file since 4636 was 2013, checked in by mervart, 15 years ago

* empty log message *
File size: 18.8 KB
Line 
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
20namespace 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
36void 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
47void 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
58void 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
69void 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
80Real 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
94void 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
166void 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
239void 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
255void 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
271void 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
287void 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
314void 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
330void 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
352void 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
419void 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
437void 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
452void 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
460void 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
470void 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
486void 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
502void 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
519void 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
536void 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
544void 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
552void 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
560void MatrixRowCol::Copy(Real r)
561{
562 // THIS = r
563 REPORT Real* elx = data; int l = storage; while (l--) *elx++ = r;
564}
565
566void MatrixRowCol::Zero()
567{
568 // THIS = 0
569 REPORT Real* elx = data; int l = storage; while (l--) *elx++ = 0;
570}
571
572void MatrixRowCol::Multiply(Real r)
573{
574 // THIS *= r
575 REPORT Real* elx = data; int l = storage; while (l--) *elx++ *= r;
576}
577
578void MatrixRowCol::Add(Real r)
579{
580 // THIS += r
581 REPORT
582 Real* elx = data; int l = storage; while (l--) *elx++ += r;
583}
584
585Real 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.
596Real 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
606Real 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
616Real 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
626Real 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
635Real 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
643void 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///@}
Note: See TracBrowser for help on using the repository browser.