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