00001
00002
00003
00004
00005
00006
00011 #include "fvar.hpp"
00012
00017 double sum(const d3_array& m)
00018 {
00019 double tmp=0.;
00020 for (int i=m.indexmin();i<=m.indexmax();i++)
00021 {
00022 tmp+=sum(m.elem(i));
00023 }
00024 return tmp;
00025 }
00026
00031 d3_array d3_array::sub(int nrl,int nrh)
00032 {
00033 if (allocated(*this))
00034 {
00035 d3_array tmp(nrl,nrh);
00036 for (int i=nrl; i<=nrh; i++)
00037 {
00038 tmp[i].shallow_copy((*this)(i));
00039 }
00040 return tmp;
00041 }
00042 else
00043 {
00044 return *this;
00045 }
00046 }
00047
00052 d3_array::d3_array(const d3_array_position& tpos)
00053 {
00054 int nrl=tpos.indexmin();
00055 int nrh=tpos.indexmax();
00056 allocate(nrl,nrh);
00057 }
00058
00063 d3_array::d3_array(int nrl,int nrh)
00064 {
00065 allocate(nrl,nrh);
00066 }
00067
00072 d3_array::d3_array(int sl,int sh,int nrl,int nrh,int ncl,int nch)
00073 {
00074 allocate(sl,sh,nrl,nrh,ncl,nch);
00075 #ifndef OPT_LIB
00076 initialize();
00077 #endif
00078 }
00079
00084 d3_array::d3_array(int sl, int sh, int nrl, int nrh, const ivector& ncl,
00085 int nch)
00086 {
00087 allocate(sl,sh,nrl,nrh,ncl,nch);
00088 #ifndef OPT_LIB
00089 initialize();
00090 #endif
00091 }
00092
00097 d3_array::d3_array(int sl, int sh, int nrl, int nrh, const ivector& ncl,
00098 const ivector& nch)
00099 {
00100 allocate(sl,sh,nrl,nrh,ncl,nch);
00101 #ifndef OPT_LIB
00102 initialize();
00103 #endif
00104 }
00105
00110 void d3_array::allocate(int sl, int sh, int nrl, int nrh, const ivector& ncl,
00111 int nch)
00112 {
00113 #ifdef DIAG
00114 myheapcheck("Entering d3_array matrix(sl,sh,nrl,nrh,ncl,nch)" );
00115 #endif
00116 if (sl !=ncl.indexmin() || sh !=ncl.indexmax())
00117 {
00118 cerr << "Incompatible array bounds in "
00119 " dmatrix(int nrl,int nrh,const ivector& ncl,const ivector& nch)" << endl;
00120 ad_exit(1);
00121 }
00122 if ( (shape=new three_array_shape(sl,sh)) == 0)
00123 {
00124 cerr << " Error allocating memory in d3_array contructor" << endl;
00125 }
00126 int ss=slicesize();
00127 if ( (t = new dmatrix[ss]) == 0)
00128 {
00129 cerr << " Error allocating memory in d3_array contructor" << endl;
00130 ad_exit(21);
00131 }
00132 t -= slicemin();
00133 for (int i=sl; i<=sh; i++)
00134 {
00135 t[i].allocate(nrl,nrh,ncl(i),nch);
00136 }
00137 }
00138
00143 void d3_array::allocate(int sl, int sh, int nrl, int nrh, int ncl,
00144 const ivector& nch)
00145 {
00146 #ifdef DIAG
00147 myheapcheck("Entering d3_array matrix(sl,sh,nrl,nrh,ncl,nch)" );
00148 #endif
00149 if (sl !=nch.indexmin() || sh !=nch.indexmax())
00150 {
00151 cerr << "Incompatible array bounds in "
00152 " dmatrix(int nrl,int nrh, const ivector& ncl, const ivector& nch)"
00153 << endl;
00154 ad_exit(1);
00155 }
00156 if ( (shape=new three_array_shape(sl,sh)) == 0)
00157 {
00158 cerr << " Error allocating memory in d3_array contructor" << endl;
00159 }
00160 int ss=slicesize();
00161 if ( (t = new dmatrix[ss]) == 0)
00162 {
00163 cerr << " Error allocating memory in d3_array contructor" << endl;
00164 ad_exit(21);
00165 }
00166 t -= slicemin();
00167 for (int i=sl; i<=sh; i++)
00168 {
00169 t[i].allocate(nrl,nrh,ncl,nch(i));
00170 }
00171 }
00172
00177 void d3_array::allocate(int sl,int sh,int nrl,int nrh,int ncl,int nch)
00178 {
00179 if ( (shape=new three_array_shape(sl,sh)) == 0)
00180 {
00181 cerr << " Error allocating memory in d3_array contructor" << endl;
00182 }
00183 int ss=slicesize();
00184 if ( (t = new dmatrix[ss]) == 0)
00185 {
00186 cerr << " Error allocating memory in d3_array contructor" << endl;
00187 ad_exit(21);
00188 }
00189 t -= slicemin();
00190 for (int i=sl; i<=sh; i++)
00191 {
00192 t[i].allocate(nrl,nrh,ncl,nch);
00193 }
00194 }
00195
00200 void d3_array::allocate(int sl,int sh,int nrl,int nrh)
00201 {
00202 if ( (shape=new three_array_shape(sl,sh)) == 0)
00203 {
00204 cerr << " Error allocating memory in d3_array contructor" << endl;
00205 }
00206 int ss=slicesize();
00207 if ( (t = new dmatrix[ss]) == 0)
00208 {
00209 cerr << " Error allocating memory in d3_array contructor" << endl;
00210 ad_exit(21);
00211 }
00212 t -= slicemin();
00213 for (int i=sl; i<=sh; i++)
00214 {
00215 t[i].allocate(nrl,nrh);
00216 }
00217 }
00218
00223 void d3_array::allocate(int sl,int sh,
00224 const index_type& nrl,const index_type& nrh)
00225 {
00226 if ( (shape=new three_array_shape(sl,sh)) == 0)
00227 {
00228 cerr << " Error allocating memory in d3_array contructor" << endl;
00229 }
00230 int ss=slicesize();
00231 if ( (t = new dmatrix[ss]) == 0)
00232 {
00233 cerr << " Error allocating memory in d3_array contructor" << endl;
00234 ad_exit(21);
00235 }
00236 t -= slicemin();
00237 for (int i=sl; i<=sh; i++)
00238 {
00239 t[i].allocate(nrl(i),nrh(i));
00240 }
00241 }
00242
00247 void d3_array::allocate(int sl,int sh)
00248 {
00249 if ( (shape=new three_array_shape(sl,sh)) == 0)
00250 {
00251 cerr << " Error allocating memory in d3_array contructor" << endl;
00252 }
00253 int ss=slicesize();
00254 if ( (t = new dmatrix[ss]) == 0)
00255 {
00256 cerr << " Error allocating memory in d3_array contructor" << endl;
00257 ad_exit(21);
00258 }
00259 t -= slicemin();
00260 for (int i=sl; i<=sh; i++)
00261 {
00262 t[i].allocate();
00263 }
00264 }
00265
00270 void d3_array::allocate(const d3_array& d3v)
00271 {
00272 int sl=d3v.slicemin();
00273 int sh=d3v.slicemax();
00274
00275
00276
00277
00278 if ( (shape=new three_array_shape(sl,sh)) == 0)
00279 {
00280 cerr << " Error allocating memory in d3_array contructor" << endl;
00281 }
00282 int ss=slicesize();
00283 if ( (t = new dmatrix[ss]) == 0)
00284 {
00285 cerr << " Error allocating memory in d3_array contructor" << endl;
00286 ad_exit(21);
00287 }
00288 t -= slicemin();
00289 for (int i=sl; i<=sh; i++)
00290 {
00291 t[i].allocate(d3v(i));
00292 }
00293 }
00294
00299 void d3_array::allocate(int sl, int sh, const ivector& nrl, const ivector& nrh,
00300 const imatrix& ncl, const imatrix& nch)
00301 {
00302 if (sl !=nrl.indexmin() || sh !=nrl.indexmax() ||
00303 sl !=nrh.indexmin() || sh !=nrh.indexmax())
00304 {
00305 cerr << "Incompatible array bounds in "
00306 " dmatrix(int nrl,int nrh,const ivector& ncl,const ivector& nch)" << endl;
00307 ad_exit(1);
00308 }
00309 if ( (shape=new three_array_shape(sl,sh)) == 0)
00310 {
00311 cerr << " Error allocating memory in d3_array contructor" << endl;
00312 }
00313 int ss=slicesize();
00314 if ( (t = new dmatrix[ss]) == 0)
00315 {
00316 cerr << " Error allocating memory in d3_array contructor" << endl;
00317 ad_exit(21);
00318 }
00319 t -= slicemin();
00320 for (int i=sl; i<=sh; i++)
00321 {
00322 t[i].allocate(nrl(i),nrh(i),ncl(i),nch(i));
00323 }
00324 }
00325
00330 void d3_array::allocate(int sl, int sh, const ivector& nrl, const ivector& nrh,
00331 int ncl, int nch)
00332 {
00333 if (sl !=nrl.indexmin() || sh !=nrl.indexmax() ||
00334 sl !=nrh.indexmin() || sh !=nrh.indexmax())
00335 {
00336 cerr << "Incompatible array bounds in "
00337 " dmatrix(int nrl,int nrh,const ivector& ncl, const ivector& nch)" << endl;
00338 ad_exit(1);
00339 }
00340 if ( (shape=new three_array_shape(sl,sh)) == 0)
00341 {
00342 cerr << " Error allocating memory in d3_array contructor" << endl;
00343 }
00344 int ss=slicesize();
00345 if ( (t = new dmatrix[ss]) == 0)
00346 {
00347 cerr << " Error allocating memory in d3_array contructor" << endl;
00348 ad_exit(21);
00349 }
00350 t -= slicemin();
00351 for (int i=sl; i<=sh; i++)
00352 {
00353 t[i].allocate(nrl(i),nrh(i),ncl,nch);
00354 }
00355 }
00356
00361 void d3_array::allocate(int sl, int sh, const ivector& nrl, int nrh, int ncl,
00362 int nch)
00363 {
00364 if (sl !=nrl.indexmin() || sh !=nrl.indexmax())
00365 {
00366 cerr << "Incompatible array bounds in "
00367 "dmatrix(int nrl,int nrh, const ivector& ncl, const ivector& nch)" << endl;
00368 ad_exit(1);
00369 }
00370 if ( (shape=new three_array_shape(sl,sh)) == 0)
00371 {
00372 cerr << " Error allocating memory in d3_array contructor" << endl;
00373 }
00374 int ss=slicesize();
00375 if ( (t = new dmatrix[ss]) == 0)
00376 {
00377 cerr << " Error allocating memory in d3_array contructor" << endl;
00378 ad_exit(21);
00379 }
00380 t -= slicemin();
00381 for (int i=sl; i<=sh; i++)
00382 {
00383 t[i].allocate(nrl(i),nrh,ncl,nch);
00384 }
00385 }
00386
00391 void d3_array::allocate(int sl, int sh, int nrl, const ivector& nrh, int ncl,
00392 int nch)
00393 {
00394 if ( sl !=nrh.indexmin() || sh !=nrh.indexmax())
00395 {
00396 cerr << "Incompatible array bounds in "
00397 "dmatrix(int nrl,int nrh,const ivector& ncl,const ivector& nch)" << endl;
00398 ad_exit(1);
00399 }
00400 if ( (shape=new three_array_shape(sl,sh)) == 0)
00401 {
00402 cerr << " Error allocating memory in d3_array contructor" << endl;
00403 }
00404 int ss=slicesize();
00405 if ( (t = new dmatrix[ss]) == 0)
00406 {
00407 cerr << " Error allocating memory in d3_array contructor" << endl;
00408 ad_exit(21);
00409 }
00410 t -= slicemin();
00411 for (int i=sl; i<=sh; i++)
00412 {
00413 t[i].allocate(nrl,nrh(i),ncl,nch);
00414 }
00415 }
00416
00421 void d3_array::allocate(int sl, int sh, const ivector& nrl, const ivector& nrh,
00422 int ncl, const imatrix& nch)
00423 {
00424 if (sl !=nrl.indexmin() || sh !=nrl.indexmax() ||
00425 sl !=nrh.indexmin() || sh !=nrh.indexmax())
00426 {
00427 cerr << "Incompatible array bounds in "
00428 "dmatrix(int nrl,int nrh,int ncl,const ivector& nch)" << endl;
00429 ad_exit(1);
00430 }
00431 if ( (shape=new three_array_shape(sl,sh)) == 0)
00432 {
00433 cerr << " Error allocating memory in d3_array contructor" << endl;
00434 }
00435 int ss=slicesize();
00436 if ( (t = new dmatrix[ss]) == 0)
00437 {
00438 cerr << " Error allocating memory in d3_array contructor" << endl;
00439 ad_exit(21);
00440 }
00441 t -= slicemin();
00442 for (int i=sl; i<=sh; i++)
00443 {
00444 t[i].allocate(nrl(i),nrh(i),ncl,nch(i));
00445 }
00446 }
00447
00452 void d3_array::allocate(int sl, int sh, int nrl, const ivector& nrh, int ncl,
00453 const imatrix& nch)
00454 {
00455 if (sl !=nrh.indexmin() || sh !=nrh.indexmax())
00456 {
00457 cerr << "Incompatible array bounds in "
00458 "dmatrix(int nrl,int nrh,int ncl,const ivector& nch)" << endl;
00459 ad_exit(1);
00460 }
00461 if ( (shape=new three_array_shape(sl,sh)) == 0)
00462 {
00463 cerr << " Error allocating memory in d3_array contructor" << endl;
00464 }
00465 int ss=slicesize();
00466 if ( (t = new dmatrix[ss]) == 0)
00467 {
00468 cerr << " Error allocating memory in d3_array contructor" << endl;
00469 ad_exit(21);
00470 }
00471 t -= slicemin();
00472 for (int i=sl; i<=sh; i++)
00473 {
00474 t[i].allocate(nrl,nrh(i),ncl,nch(i));
00475 }
00476 }
00477
00482 void d3_array::allocate(int sl, int sh, const ivector& nrl, const ivector& nrh,
00483 const ivector& ncl, const ivector& nch)
00484 {
00485 if (sl !=nrl.indexmin() || sh !=nrl.indexmax() ||
00486 sl !=nrh.indexmin() || sh !=nrh.indexmax())
00487 {
00488 cerr << "Incompatible array bounds in "
00489 "dmatrix(int nrl,int nrh,const ivector& ncl,const ivector& nch)" << endl;
00490 ad_exit(1);
00491 }
00492 if ( (shape=new three_array_shape(sl,sh)) == 0)
00493 {
00494 cerr << " Error allocating memory in d3_array contructor" << endl;
00495 }
00496 int ss=slicesize();
00497 if ( (t = new dmatrix[ss]) == 0)
00498 {
00499 cerr << " Error allocating memory in d3_array contructor" << endl;
00500 ad_exit(21);
00501 }
00502 t -= slicemin();
00503 for (int i=sl; i<=sh; i++)
00504 {
00505 t[i].allocate(nrl(i),nrh(i),ncl(i),nch(i));
00506 }
00507 }
00508
00513 void d3_array::allocate(int sl, int sh, int nrl, int nrh, const ivector& ncl,
00514 const ivector& nch)
00515 {
00516 if ( (shape=new three_array_shape(sl,sh)) == 0)
00517 {
00518 cerr << " Error allocating memory in d3_array contructor" << endl;
00519 }
00520 int ss=slicesize();
00521 if ( (t = new dmatrix[ss]) == 0)
00522 {
00523 cerr << " Error allocating memory in d3_array contructor" << endl;
00524 ad_exit(21);
00525 }
00526 t -= slicemin();
00527 for (int i=sl; i<=sh; i++)
00528 {
00529 t[i].allocate(nrl,nrh,ncl(i),nch(i));
00530 }
00531 }
00532
00537 void d3_array::allocate(int sl, int sh, int nrl, const ivector& nrh, int ncl,
00538 const ivector& nch)
00539 {
00540 if ( (shape=new three_array_shape(sl,sh)) == 0)
00541 {
00542 cerr << " Error allocating memory in d3_array contructor" << endl;
00543 }
00544 int ss=slicesize();
00545 if ( (t = new dmatrix[ss]) == 0)
00546 {
00547 cerr << " Error allocating memory in d3_array contructor" << endl;
00548 ad_exit(21);
00549 }
00550 t -= slicemin();
00551 for (int i=sl; i<=sh; i++)
00552 {
00553 t[i].allocate(nrl,nrh(i),ncl,nch(i));
00554 }
00555 }
00556
00561 d3_array::d3_array(int sl, int sh, int nrl, const ivector& nrh, int ncl,
00562 const ivector& nch)
00563 {
00564 allocate(sl,sh,nrl,nrh,ncl,nch);
00565 }
00566
00571 d3_array::d3_array(int sl, int sh, const ivector& nrl, const ivector& nrh,
00572 const imatrix& ncl, const imatrix& nch)
00573 {
00574 allocate(sl,sh,nrl,nrh,ncl,nch);
00575 }
00576
00581 d3_array::d3_array(int sl, int sh, const ivector& nrl, const ivector& nrh,
00582 int ncl, const imatrix& nch)
00583 {
00584 allocate(sl,sh,nrl,nrh,ncl,nch);
00585 }
00586
00591 d3_array::d3_array(int sl, int sh, int nrl, const ivector& nrh, int ncl,
00592 const imatrix& nch)
00593 {
00594 allocate(sl,sh,nrl,nrh,ncl,nch);
00595 }
00596
00601 d3_array::d3_array(int sl, int sh, int nrl, const ivector& nrh, int ncl,
00602 int nch)
00603 {
00604 allocate(sl,sh,nrl,nrh,ncl,nch);
00605 }
00606
00611 d3_array::d3_array(const d3_array& m2)
00612 {
00613 shape=m2.shape;
00614 if (shape)
00615 {
00616 (shape->ncopies)++;
00617 }
00618 #ifdef SAFE_ALL
00619 else
00620 {
00621 cerr << "Making a copy of an unallocated d3_array"<<endl;
00622 }
00623 #endif
00624 t = m2.t;
00625 }
00626
00631 void d3_array::shallow_copy(const d3_array& m2)
00632 {
00633 shape=m2.shape;
00634 if (shape)
00635 {
00636 (shape->ncopies)++;
00637 }
00638 #ifdef SAFE_ALL
00639 else
00640 {
00641 cerr << "Making a copy of an unallocated d3_array"<<endl;
00642 }
00643 #endif
00644 t = m2.t;
00645 }
00646
00651 void d3_array::initialize()
00652 {
00653 if (!(!(*this)))
00654 {
00655 for (int i=slicemin();i<=slicemax();i++)
00656 {
00657 elem(i).initialize();
00658 }
00659 }
00660 }
00661
00666 void d3_array::deallocate()
00667 {
00668 if (shape)
00669 {
00670
00671
00672 t += slicemin();
00673
00674 delete [] t;
00675
00676 delete shape;
00677
00678 t=NULL;
00679 shape=NULL;
00680 }
00681 #ifdef SAFE_ALL
00682 else
00683 {
00684 cerr << "Warning -- trying to deallocate an unallocated dmatrix"<<endl;
00685 }
00686 #endif
00687 }
00688
00692 d3_array::~d3_array()
00693 {
00694 if (shape)
00695 {
00696 if (shape->ncopies)
00697 {
00698 (shape->ncopies)--;
00699 }
00700 else
00701 {
00702 deallocate();
00703 }
00704 }
00705 #ifdef SAFE_ALL
00706 else
00707 {
00708 cerr << "Warning -- trying to deallocate an unallocated d3_array"<<endl;
00709 }
00710 #endif
00711 }
00716 three_array_shape::three_array_shape(int sl,int su)
00717 {
00718 slice_min=sl;
00719 slice_max=su;
00720
00721
00722
00723
00724
00725
00726
00727 ncopies=0;
00728 }