Go to the documentation of this file.00001
00002
00003
00008
00009
00010
00011
00012 #define WANT_MATH
00013
00014 #include "include.h"
00015
00016 #include "newmat.h"
00017 #include "newmatrc.h"
00018
00019 #ifdef use_namespace
00020 namespace NEWMAT {
00021 #endif
00022
00023
00024 #ifdef DO_REPORT
00025 #define REPORT { static ExeCounter ExeCount(__LINE__,2); ++ExeCount; }
00026 #else
00027 #define REPORT {}
00028 #endif
00029
00030
00031
00032 #define MONITOR(what,store,storage) {}
00033
00034
00035
00036 void MatrixRowCol::Add(const MatrixRowCol& mrc)
00037 {
00038
00039 REPORT
00040 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00041 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00042 if (l<=0) return;
00043 Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00044 while (l--) *elx++ += *el++;
00045 }
00046
00047 void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x)
00048 {
00049 REPORT
00050
00051 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00052 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00053 if (l<=0) return;
00054 Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00055 while (l--) *elx++ += *el++ * x;
00056 }
00057
00058 void MatrixRowCol::Sub(const MatrixRowCol& mrc)
00059 {
00060 REPORT
00061
00062 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00063 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00064 if (l<=0) return;
00065 Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00066 while (l--) *elx++ -= *el++;
00067 }
00068
00069 void MatrixRowCol::Inject(const MatrixRowCol& mrc)
00070
00071 {
00072 REPORT
00073 int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00074 if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00075 if (l<=0) return;
00076 Real* elx=data+(f-skip); Real* ely=mrc.data+(f-mrc.skip);
00077 while (l--) *elx++ = *ely++;
00078 }
00079
00080 Real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00081 {
00082 REPORT
00083 int f = mrc1.skip; int f2 = mrc2.skip;
00084 int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
00085 if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
00086 if (l<=0) return 0.0;
00087
00088 Real* el1=mrc1.data+(f-mrc1.skip); Real* el2=mrc2.data+(f-mrc2.skip);
00089 Real sum = 0.0;
00090 while (l--) sum += *el1++ * *el2++;
00091 return sum;
00092 }
00093
00094 void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00095 {
00096
00097 int f = skip; int l = skip + storage;
00098 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00099 if (f1<f) f1=f; if (l1>l) l1=l;
00100 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00101 if (f2<f) f2=f; if (l2>l) l2=l;
00102 Real* el = data + (f-skip);
00103 Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
00104 if (f1<f2)
00105 {
00106 int i = f1-f; while (i--) *el++ = 0.0;
00107 if (l1<=f2)
00108 {
00109 REPORT
00110 i = l1-f1; while (i--) *el++ = *el1++;
00111 i = f2-l1; while (i--) *el++ = 0.0;
00112 i = l2-f2; while (i--) *el++ = *el2++;
00113 i = l-l2; while (i--) *el++ = 0.0;
00114 }
00115 else
00116 {
00117 i = f2-f1; while (i--) *el++ = *el1++;
00118 if (l1<=l2)
00119 {
00120 REPORT
00121 i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
00122 i = l2-l1; while (i--) *el++ = *el2++;
00123 i = l-l2; while (i--) *el++ = 0.0;
00124 }
00125 else
00126 {
00127 REPORT
00128 i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
00129 i = l1-l2; while (i--) *el++ = *el1++;
00130 i = l-l1; while (i--) *el++ = 0.0;
00131 }
00132 }
00133 }
00134 else
00135 {
00136 int i = f2-f; while (i--) *el++ = 0.0;
00137 if (l2<=f1)
00138 {
00139 REPORT
00140 i = l2-f2; while (i--) *el++ = *el2++;
00141 i = f1-l2; while (i--) *el++ = 0.0;
00142 i = l1-f1; while (i--) *el++ = *el1++;
00143 i = l-l1; while (i--) *el++ = 0.0;
00144 }
00145 else
00146 {
00147 i = f1-f2; while (i--) *el++ = *el2++;
00148 if (l2<=l1)
00149 {
00150 REPORT
00151 i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
00152 i = l1-l2; while (i--) *el++ = *el1++;
00153 i = l-l1; while (i--) *el++ = 0.0;
00154 }
00155 else
00156 {
00157 REPORT
00158 i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
00159 i = l2-l1; while (i--) *el++ = *el2++;
00160 i = l-l2; while (i--) *el++ = 0.0;
00161 }
00162 }
00163 }
00164 }
00165
00166 void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00167 {
00168
00169 int f = skip; int l = skip + storage;
00170 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00171 if (f1<f) f1=f; if (l1>l) l1=l;
00172 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00173 if (f2<f) f2=f; if (l2>l) l2=l;
00174 Real* el = data + (f-skip);
00175 Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
00176 if (f1<f2)
00177 {
00178 int i = f1-f; while (i--) *el++ = 0.0;
00179 if (l1<=f2)
00180 {
00181 REPORT
00182 i = l1-f1; while (i--) *el++ = *el1++;
00183 i = f2-l1; while (i--) *el++ = 0.0;
00184 i = l2-f2; while (i--) *el++ = - *el2++;
00185 i = l-l2; while (i--) *el++ = 0.0;
00186 }
00187 else
00188 {
00189 i = f2-f1; while (i--) *el++ = *el1++;
00190 if (l1<=l2)
00191 {
00192 REPORT
00193 i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
00194 i = l2-l1; while (i--) *el++ = - *el2++;
00195 i = l-l2; while (i--) *el++ = 0.0;
00196 }
00197 else
00198 {
00199 REPORT
00200 i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
00201 i = l1-l2; while (i--) *el++ = *el1++;
00202 i = l-l1; while (i--) *el++ = 0.0;
00203 }
00204 }
00205 }
00206 else
00207 {
00208 int i = f2-f; while (i--) *el++ = 0.0;
00209 if (l2<=f1)
00210 {
00211 REPORT
00212 i = l2-f2; while (i--) *el++ = - *el2++;
00213 i = f1-l2; while (i--) *el++ = 0.0;
00214 i = l1-f1; while (i--) *el++ = *el1++;
00215 i = l-l1; while (i--) *el++ = 0.0;
00216 }
00217 else
00218 {
00219 i = f1-f2; while (i--) *el++ = - *el2++;
00220 if (l2<=l1)
00221 {
00222 REPORT
00223 i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
00224 i = l1-l2; while (i--) *el++ = *el1++;
00225 i = l-l1; while (i--) *el++ = 0.0;
00226 }
00227 else
00228 {
00229 REPORT
00230 i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
00231 i = l2-l1; while (i--) *el++ = - *el2++;
00232 i = l-l2; while (i--) *el++ = 0.0;
00233 }
00234 }
00235 }
00236 }
00237
00238
00239 void MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x)
00240 {
00241
00242 REPORT
00243 if (!storage) return;
00244 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00245 if (f < skip) { f = skip; if (l < f) l = f; }
00246 if (l > lx) { l = lx; if (f > lx) f = lx; }
00247
00248 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00249
00250 int l1 = f-skip; while (l1--) *elx++ = x;
00251 l1 = l-f; while (l1--) *elx++ = *ely++ + x;
00252 lx -= l; while (lx--) *elx++ = x;
00253 }
00254
00255 void MatrixRowCol::NegAdd(const MatrixRowCol& mrc1, Real x)
00256 {
00257
00258 REPORT
00259 if (!storage) return;
00260 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00261 if (f < skip) { f = skip; if (l < f) l = f; }
00262 if (l > lx) { l = lx; if (f > lx) f = lx; }
00263
00264 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00265
00266 int l1 = f-skip; while (l1--) *elx++ = x;
00267 l1 = l-f; while (l1--) *elx++ = x - *ely++;
00268 lx -= l; while (lx--) *elx++ = x;
00269 }
00270
00271 void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
00272 {
00273
00274 REPORT
00275 if (!storage) return;
00276 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00277 if (f < skip) { f = skip; if (l < f) l = f; }
00278 if (l > lx) { l = lx; if (f > lx) f = lx; }
00279
00280 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00281
00282 int l1 = f-skip; while (l1--) { *elx = - *elx; elx++; }
00283 l1 = l-f; while (l1--) { *elx = *ely++ - *elx; elx++; }
00284 lx -= l; while (lx--) { *elx = - *elx; elx++; }
00285 }
00286
00287 void MatrixRowCol::ConCat(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00288 {
00289
00290 REPORT
00291 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage; int lx = skip + storage;
00292 if (f1 < skip) { f1 = skip; if (l1 < f1) l1 = f1; }
00293 if (l1 > lx) { l1 = lx; if (f1 > lx) f1 = lx; }
00294
00295 Real* elx = data;
00296
00297 int i = f1-skip; while (i--) *elx++ =0.0;
00298 i = l1-f1;
00299 if (i)
00300 { Real* ely = mrc1.data+(f1-mrc1.skip); while (i--) *elx++ = *ely++; }
00301
00302 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage; i = mrc1.length;
00303 int skipx = l1 - i; lx -= i;
00304 if (f2 < skipx) { f2 = skipx; if (l2 < f2) l2 = f2; }
00305 if (l2 > lx) { l2 = lx; if (f2 > lx) f2 = lx; }
00306
00307 i = f2-skipx; while (i--) *elx++ = 0.0;
00308 i = l2-f2;
00309 if (i)
00310 { Real* ely = mrc2.data+(f2-mrc2.skip); while (i--) *elx++ = *ely++; }
00311 lx -= l2; while (lx--) *elx++ = 0.0;
00312 }
00313
00314 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1)
00315
00316 {
00317 REPORT
00318 if (!storage) return;
00319 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00320 if (f < skip) { f = skip; if (l < f) l = f; }
00321 if (l > lx) { l = lx; if (f > lx) f = lx; }
00322
00323 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00324
00325 int l1 = f-skip; while (l1--) *elx++ = 0;
00326 l1 = l-f; while (l1--) *elx++ *= *ely++;
00327 lx -= l; while (lx--) *elx++ = 0;
00328 }
00329
00330 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00331
00332 {
00333 int f = skip; int l = skip + storage;
00334 int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00335 if (f1<f) f1=f; if (l1>l) l1=l;
00336 int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00337 if (f2<f) f2=f; if (l2>l) l2=l;
00338 Real* el = data + (f-skip); int i;
00339 if (f1<f2) f1 = f2; if (l1>l2) l1 = l2;
00340 if (l1<=f1) { REPORT i = l-f; while (i--) *el++ = 0.0; }
00341 else
00342 {
00343 REPORT
00344 Real* el1 = mrc1.data+(f1-mrc1.skip);
00345 Real* el2 = mrc2.data+(f1-mrc2.skip);
00346 i = f1-f ; while (i--) *el++ = 0.0;
00347 i = l1-f1; while (i--) *el++ = *el1++ * *el2++;
00348 i = l-l1; while (i--) *el++ = 0.0;
00349 }
00350 }
00351
00352 void MatrixRowCol::KP(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00353
00354 {
00355 int f = skip; int s = storage; Real* el = data; int i;
00356
00357 i = mrc1.skip * mrc2.length;
00358 if (i > f)
00359 {
00360 i -= f; f = 0; if (i > s) { i = s; s = 0; } else s -= i;
00361 while (i--) *el++ = 0.0;
00362 if (s == 0) return;
00363 }
00364 else f -= i;
00365
00366 i = mrc1.storage; Real* el1 = mrc1.data;
00367 int mrc2_skip = mrc2.skip; int mrc2_storage = mrc2.storage;
00368 int mrc2_length = mrc2.length;
00369 int mrc2_remain = mrc2_length - mrc2_skip - mrc2_storage;
00370 while (i--)
00371 {
00372 int j; Real* el2 = mrc2.data; Real vel1 = *el1;
00373 if (f == 0 && mrc2_length <= s)
00374 {
00375 j = mrc2_skip; s -= j; while (j--) *el++ = 0.0;
00376 j = mrc2_storage; s -= j; while (j--) *el++ = vel1 * *el2++;
00377 j = mrc2_remain; s -= j; while (j--) *el++ = 0.0;
00378 }
00379 else if (f >= mrc2_length) f -= mrc2_length;
00380 else
00381 {
00382 j = mrc2_skip;
00383 if (j > f)
00384 {
00385 j -= f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
00386 while (j--) *el++ = 0.0;
00387 }
00388 else f -= j;
00389
00390 j = mrc2_storage;
00391 if (j > f)
00392 {
00393 j -= f; el2 += f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
00394 while (j--) *el++ = vel1 * *el2++;
00395 }
00396 else f -= j;
00397
00398 j = mrc2_remain;
00399 if (j > f)
00400 {
00401 j -= f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
00402 while (j--) *el++ = 0.0;
00403 }
00404 else f -= j;
00405 }
00406 if (s == 0) return;
00407 ++el1;
00408 }
00409
00410 i = (mrc1.length - mrc1.skip - mrc1.storage) * mrc2.length;
00411 if (i > f)
00412 {
00413 i -= f; if (i > s) i = s;
00414 while (i--) *el++ = 0.0;
00415 }
00416 }
00417
00418
00419 void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
00420 {
00421
00422 REPORT
00423 if (!storage) return;
00424 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00425 if (f < skip) { f = skip; if (l < f) l = f; }
00426 if (l > lx) { l = lx; if (f > lx) f = lx; }
00427
00428 Real* elx = data; Real* ely = 0;
00429
00430 if (l-f) ely = mrc1.data+(f-mrc1.skip);
00431
00432 int l1 = f-skip; while (l1--) *elx++ = 0.0;
00433 l1 = l-f; while (l1--) *elx++ = *ely++;
00434 lx -= l; while (lx--) *elx++ = 0.0;
00435 }
00436
00437 void MatrixRowCol::CopyCheck(const MatrixRowCol& mrc1)
00438
00439 {
00440 REPORT
00441 if (!storage) return;
00442 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00443 if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
00444
00445 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00446
00447 int l1 = f-skip; while (l1--) *elx++ = 0.0;
00448 l1 = l-f; while (l1--) *elx++ = *ely++;
00449 lx -= l; while (lx--) *elx++ = 0.0;
00450 }
00451
00452 void MatrixRowCol::Check(const MatrixRowCol& mrc1)
00453
00454 {
00455 REPORT
00456 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00457 if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
00458 }
00459
00460 void MatrixRowCol::Check()
00461
00462
00463
00464 {
00465 REPORT
00466 if (skip!=0 || storage!=length)
00467 Throw(ProgramException("Illegal Conversion"));
00468 }
00469
00470 void MatrixRowCol::Negate(const MatrixRowCol& mrc1)
00471 {
00472
00473 REPORT
00474 if (!storage) return;
00475 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00476 if (f < skip) { f = skip; if (l < f) l = f; }
00477 if (l > lx) { l = lx; if (f > lx) f = lx; }
00478
00479 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00480
00481 int l1 = f-skip; while (l1--) *elx++ = 0.0;
00482 l1 = l-f; while (l1--) *elx++ = - *ely++;
00483 lx -= l; while (lx--) *elx++ = 0.0;
00484 }
00485
00486 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s)
00487 {
00488
00489 REPORT
00490 if (!storage) return;
00491 int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00492 if (f < skip) { f = skip; if (l < f) l = f; }
00493 if (l > lx) { l = lx; if (f > lx) f = lx; }
00494
00495 Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00496
00497 int l1 = f-skip; while (l1--) *elx++ = 0.0;
00498 l1 = l-f; while (l1--) *elx++ = *ely++ * s;
00499 lx -= l; while (lx--) *elx++ = 0.0;
00500 }
00501
00502 void DiagonalMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1)
00503 {
00504
00505 REPORT
00506 int f = mrc1.skip; int f0 = mrc.skip;
00507 int l = f + mrc1.storage; int lx = f0 + mrc.storage;
00508 if (f < f0) { f = f0; if (l < f) l = f; }
00509 if (l > lx) { l = lx; if (f > lx) f = lx; }
00510
00511 Real* elx = mrc.data; Real* eld = store+f;
00512
00513 int l1 = f-f0; while (l1--) *elx++ = 0.0;
00514 l1 = l-f; while (l1--) *elx++ /= *eld++;
00515 lx -= l; while (lx--) *elx++ = 0.0;
00516
00517 }
00518
00519 void IdentityMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1)
00520 {
00521
00522 REPORT
00523 int f = mrc1.skip; int f0 = mrc.skip;
00524 int l = f + mrc1.storage; int lx = f0 + mrc.storage;
00525 if (f < f0) { f = f0; if (l < f) l = f; }
00526 if (l > lx) { l = lx; if (f > lx) f = lx; }
00527
00528 Real* elx = mrc.data; Real eldv = *store;
00529
00530 int l1 = f-f0; while (l1--) *elx++ = 0.0;
00531 l1 = l-f; while (l1--) *elx++ /= eldv;
00532 lx -= l; while (lx--) *elx++ = 0.0;
00533
00534 }
00535
00536 void MatrixRowCol::Copy(const double*& r)
00537 {
00538
00539 REPORT
00540 Real* elx = data; const double* ely = r+skip; r += length;
00541 int l = storage; while (l--) *elx++ = (Real)*ely++;
00542 }
00543
00544 void MatrixRowCol::Copy(const float*& r)
00545 {
00546
00547 REPORT
00548 Real* elx = data; const float* ely = r+skip; r += length;
00549 int l = storage; while (l--) *elx++ = (Real)*ely++;
00550 }
00551
00552 void MatrixRowCol::Copy(const int*& r)
00553 {
00554
00555 REPORT
00556 Real* elx = data; const int* ely = r+skip; r += length;
00557 int l = storage; while (l--) *elx++ = (Real)*ely++;
00558 }
00559
00560 void MatrixRowCol::Copy(Real r)
00561 {
00562
00563 REPORT Real* elx = data; int l = storage; while (l--) *elx++ = r;
00564 }
00565
00566 void MatrixRowCol::Zero()
00567 {
00568
00569 REPORT Real* elx = data; int l = storage; while (l--) *elx++ = 0;
00570 }
00571
00572 void MatrixRowCol::Multiply(Real r)
00573 {
00574
00575 REPORT Real* elx = data; int l = storage; while (l--) *elx++ *= r;
00576 }
00577
00578 void MatrixRowCol::Add(Real r)
00579 {
00580
00581 REPORT
00582 Real* elx = data; int l = storage; while (l--) *elx++ += r;
00583 }
00584
00585 Real MatrixRowCol::SumAbsoluteValue()
00586 {
00587 REPORT
00588 Real sum = 0.0; Real* elx = data; int l = storage;
00589 while (l--) sum += fabs(*elx++);
00590 return sum;
00591 }
00592
00593
00594
00595
00596 Real MatrixRowCol::MaximumAbsoluteValue1(Real r, int& i)
00597 {
00598 REPORT
00599 Real* elx = data; int l = storage; int li = -1;
00600 while (l--) { Real f = fabs(*elx++); if (r <= f) { r = f; li = l; } }
00601 i = (li >= 0) ? storage - li + skip : 0;
00602 return r;
00603 }
00604
00605
00606 Real MatrixRowCol::MinimumAbsoluteValue1(Real r, int& i)
00607 {
00608 REPORT
00609 Real* elx = data; int l = storage; int li = -1;
00610 while (l--) { Real f = fabs(*elx++); if (r >= f) { r = f; li = l; } }
00611 i = (li >= 0) ? storage - li + skip : 0;
00612 return r;
00613 }
00614
00615
00616 Real MatrixRowCol::Maximum1(Real r, int& i)
00617 {
00618 REPORT
00619 Real* elx = data; int l = storage; int li = -1;
00620 while (l--) { Real f = *elx++; if (r <= f) { r = f; li = l; } }
00621 i = (li >= 0) ? storage - li + skip : 0;
00622 return r;
00623 }
00624
00625
00626 Real MatrixRowCol::Minimum1(Real r, int& i)
00627 {
00628 REPORT
00629 Real* elx = data; int l = storage; int li = -1;
00630 while (l--) { Real f = *elx++; if (r >= f) { r = f; li = l; } }
00631 i = (li >= 0) ? storage - li + skip : 0;
00632 return r;
00633 }
00634
00635 Real MatrixRowCol::Sum()
00636 {
00637 REPORT
00638 Real sum = 0.0; Real* elx = data; int l = storage;
00639 while (l--) sum += *elx++;
00640 return sum;
00641 }
00642
00643 void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
00644 {
00645 mrc.length = l1; int d = skip - skip1;
00646 if (d<0) { mrc.skip = 0; mrc.data = data - d; }
00647 else { mrc.skip = d; mrc.data = data; }
00648 d = skip + storage - skip1;
00649 d = ((l1 < d) ? l1 : d) - mrc.skip; mrc.storage = (d < 0) ? 0 : d;
00650 mrc.cw = 0;
00651 }
00652
00653 #ifdef use_namespace
00654 }
00655 #endif
00656
00657