nDarray.cppGo to the documentation of this file.00001 // $Revision: 1.6 $ 00002 // $Date: 2006/08/07 22:17:27 $ 00003 // $Source: /usr/local/cvs/OpenSees/SRC/nDarray/nDarray.cpp,v $ 00004 00005 00006 //############################################################################ 00007 //# # 00008 //# /~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/~~\ # 00009 //# | |____| # 00010 //# | | # 00011 //# | | # 00012 //# | B A S E | # 00013 //# | | # 00014 //# | | # 00015 //# | C L A S S E S | # 00016 //# | | # 00017 //# | | # 00018 //# | C + + S O U R C E | # 00019 //# | | # 00020 //# | | # 00021 //# | | # 00022 //# | | # 00023 //# | | # 00024 //# /~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/ | # 00025 //# | | | # 00026 //# \_________________________________________\__/ # 00027 //# # 00028 //# # 00029 //############################################################################ 00030 // 00031 // "C makes it easy to shoot yourself in the foot, C++ makes it harder, 00032 // but when you do, it blows away your whole leg" -- Bjarne Stroustrup 00033 // 00034 //############################################################################# 00035 // PAZI 'VAMO: 00036 // 00037 // 00038 //0. napravi vise vrsta za (int), (int, int), 00039 // (int,int,int) . . . 00040 // 00041 //1. for efficiency code val(int), val(int, int) . . . up to 00042 // fourth order inline . . . also operator() . . . 00043 // Check in Ellis & Stroustrup about inline functions 00044 // 00045 //2. Code symetric operators ( *, +, -, : . . . ) as friends 00046 // see Coplien's recomendation . . . 00047 // 00048 // 00049 // 00050 // 00053 //################################################################################ 00054 //# COPY-YES (C): :-)) # 00055 //# PROJECT: Object Oriented Finite Element Program # 00056 //# PURPOSE: # 00057 //# CLASS: nDarray # 00058 //# # 00059 //# VERSION: # 00060 //# LANGUAGE: C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 ) # 00061 //# TARGET OS: DOS || UNIX || . . . # 00062 //# DESIGNER(S): Boris Jeremic # 00063 //# PROGRAMMER(S): Boris Jeremic # 00064 //# # 00065 //# # 00066 //# DATE: May 28. - July 21 '93 # 00067 //# UPDATE HISTORY: august 05. - august __ '93 base_03 ( BJmatrix, BJvector, # 00068 //# skyBJmatrix ...) from this one # 00069 //# August 22-29 '94 choped to separate files and worked on # 00070 //# const and & issues # 00071 //# August 30-31 '94 added use_def_dim to full the CC # 00072 //# resolved problem with temoraries for # 00073 //# operators + and - ( +=, -= ) # 00074 //# January 16 '95 fixed the memory leakage introduced # 00075 //# by previous work on +=, -+. I was # 00076 //# by mistake decreasing # 00077 //# this->pc_nDarray_rep->total_numb--; # 00078 //# inststead of # 00079 //# this->pc_nDarray_rep->n--; # 00080 //# 28June2004 added val4 for efficiency still # 00081 //# to be worked on # 00082 //# # 00083 //# # 00084 //################################################################################ 00085 //*/ 00086 // nDarray.cc 00087 00088 #ifndef NDARRAY_CC 00089 #define NDARRAY_CC 00090 00091 00092 // #include "basics.hh" 00093 #include "nDarray.h" 00094 00095 //############################################################################## 00096 nDarray::nDarray(int rank_of_nDarray, double initval) 00097 { 00098 // create the structure: 00099 pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 1 00100 pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00101 00102 // in the case of nDarray_rank=0 add one to get right thing from the 00103 // operator new 00104 int one_or0 = 0; 00105 if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00106 pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00107 // dimensions 00108 const int default_dim = 1; 00109 pc_nDarray_rep->total_numb = 1; 00110 for( int idim = 0 ; idim < pc_nDarray_rep->nDarray_rank ; idim++ ) 00111 { 00112 pc_nDarray_rep->dim[idim] = default_dim; 00113 pc_nDarray_rep->total_numb *= pc_nDarray_rep->dim[idim]; 00114 } 00115 00116 // allocate memory for the actual nDarray as nDarray 00117 pc_nDarray_rep->pd_nDdata = new double [(size_t) pc_nDarray_rep->total_numb]; 00118 if (!pc_nDarray_rep->pd_nDdata) 00119 { 00120 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00121 ::exit(1); 00122 } 00123 00124 pc_nDarray_rep->n = 1; // so far, there's one reference 00125 00126 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00127 pc_nDarray_rep->pd_nDdata[i] = initval; 00128 00129 } 00130 00131 //############################################################################## 00132 nDarray::nDarray(int rank_of_nDarray, const int *pdim, double *values) 00133 { 00134 // create the structure: 00135 pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 2 00136 pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00137 00138 // in the case of nDarray_rank=0 add one to get right thing from the 00139 // operator new 00140 int one_or0 = 0; 00141 if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00142 pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00143 // dimensions 00144 00145 pc_nDarray_rep->total_numb = 1; 00146 for( int idim = 0 ; idim < pc_nDarray_rep->nDarray_rank ; idim++ ) 00147 { 00148 pc_nDarray_rep->dim[idim] = pdim[idim]; 00149 pc_nDarray_rep->total_numb *= pc_nDarray_rep->dim[idim]; 00150 } 00151 00152 00153 // allocate memory for the actual nDarray as nDarray 00154 pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00155 if (!pc_nDarray_rep->pd_nDdata) 00156 { 00157 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00158 ::exit(1); 00159 } 00160 00161 00162 pc_nDarray_rep->n = 1; // so far, there's one reference 00163 00164 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00165 pc_nDarray_rep->pd_nDdata[i] = values[i]; 00166 00167 } 00168 00169 //############################################################################## 00170 nDarray::nDarray(int rank_of_nDarray, const int *pdim, double initvalue) 00171 { 00172 // create the structure: 00173 pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 3 00174 pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00175 00176 // in the case of nDarray_rank=0 add one to get right thing from the 00177 // operator new 00178 int one_or0 = 0; 00179 if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00180 pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00181 // dimensions 00182 00183 pc_nDarray_rep->total_numb = 1; 00184 for( int idim = 0 ; idim < pc_nDarray_rep->nDarray_rank ; idim++ ) 00185 { 00186 pc_nDarray_rep->dim[idim] = pdim[idim]; 00187 pc_nDarray_rep->total_numb *= pc_nDarray_rep->dim[idim]; 00188 } 00189 00190 // allocate memory for the actual nDarray as nDarray 00191 pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00192 if (!pc_nDarray_rep->pd_nDdata) 00193 { 00194 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00195 ::exit(1); 00196 } 00197 00198 pc_nDarray_rep->n = 1; // so far, there's one reference 00199 00200 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00201 pc_nDarray_rep->pd_nDdata[i] = initvalue; 00202 } 00203 00204 00205 00206 //############################################################################## 00207 // special case for BJmatrix and BJvector . . . 00208 nDarray::nDarray(int rank_of_nDarray, int rows, int cols, double *values) 00209 { 00210 // create the structure: 00211 pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 4 00212 pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00213 00214 // not needed for BJmatrix or BJvector but who knows # 00215 // in the case of nDarray_rank=0 add one to get right thing from the 00216 // operator new 00217 int one_or0 = 0; 00218 if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00219 pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00220 // dimensions 00221 00222 pc_nDarray_rep->total_numb = 1; 00223 00224 pc_nDarray_rep->dim[0] = rows; 00225 pc_nDarray_rep->dim[1] = cols; 00226 pc_nDarray_rep->total_numb = rows*cols; 00227 00228 // allocate memory for the actual nDarray as nDarray 00229 pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00230 if (!pc_nDarray_rep->pd_nDdata) 00231 { 00232 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00233 ::exit(1); 00234 } 00235 pc_nDarray_rep->n = 1; // so far, there's one reference 00236 00237 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00238 pc_nDarray_rep->pd_nDdata[i] = values[i]; 00239 } 00240 00241 //############################################################################## 00242 // special case for BJmatrix and BJvector . . . 00243 nDarray::nDarray(int rank_of_nDarray, int rows, int cols, double values) 00244 { 00245 // create the structure: 00246 pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 5 00247 pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00248 00249 // not needed for BJmatrix or BJvector but who knows # 00250 // in the case of nDarray_rank=0 add one to get right thing from the 00251 // operator new 00252 int one_or0 = 0; 00253 if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00254 pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00255 // dimensions 00256 00257 pc_nDarray_rep->total_numb = 1; 00258 00259 pc_nDarray_rep->dim[0] = rows; 00260 pc_nDarray_rep->dim[1] = cols; 00261 pc_nDarray_rep->total_numb = rows*cols; 00262 00263 // allocate memory for the actual nDarray as nDarray 00264 pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00265 if (!pc_nDarray_rep->pd_nDdata) 00266 { 00267 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00268 ::exit(1); 00269 } 00270 pc_nDarray_rep->n = 1; // so far, there's one reference 00271 00272 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00273 pc_nDarray_rep->pd_nDdata[i] = values; 00274 } 00275 00276 00277 //---//############################################################################## 00278 //---// special case for BJmatrix and BJvector . . . 00279 //---nDarray::nDarray(int rank_of_nDarray, int rows, int cols, double initvalue) 00280 //---{ 00281 //--- // create the structure: 00282 //--- pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 00283 //--- pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00284 //--- 00285 //---// not needed for BJmatrix or BJvector but who knows # 00286 //---// in the case of nDarray_rank=0 add one to get right thing from the 00287 //---// operator new 00288 //--- int one_or0 = 0; 00289 //--- if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00290 //--- pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00291 //--- // dimensions 00292 //--- 00293 //--- pc_nDarray_rep->total_numb = 1; 00294 //--- 00295 //--- pc_nDarray_rep->dim[0] = rows; 00296 //--- pc_nDarray_rep->dim[1] = cols; 00297 //--- pc_nDarray_rep->total_numb = rows*cols; 00298 //--- 00299 //---// allocate memory for the actual nDarray as nDarray 00300 //--- pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00301 //--- if (!pc_nDarray_rep->pd_nDdata) 00302 //--- { 00303 //--- ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00304 //--- ::exit(1); 00305 //--- } 00306 //--- 00307 //--- pc_nDarray_rep->n = 1; // so far, there's one reference 00308 //--- 00309 //--- for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00310 //--- pc_nDarray_rep->pd_nDdata[i] = initvalue; 00311 //---} 00312 //--- 00313 00314 00315 00316 00317 00318 //...//############################################################################## 00319 // create a unit nDarray 00320 nDarray::nDarray(const char *flag, int rank_of_nDarray, const int *pdim) 00321 { 00322 if ( flag[0] != 'I' && flag[0] != 'e' && flag[0] != 'C' ) 00323 { 00324 ::printf("\nTo create a 2nd rank Kronecker delta type: nDarray (\"I\",2,dims);\n"); 00325 // ::printf( "To create a 4th rank unit nDarray type: nDarray (\"I\",4,dims);\n"); 00326 ::printf( "To create a 3th rank Levi-Civita BJtensor type: nDarray (\"e\",3,dims);\n"); 00327 ::printf( "To create a 2nd rank Cosserat Kronecker delta type: nDarray (\"C\",3,dims);\n"); 00328 ::exit( 1 ); 00329 } 00330 // create the structure: 00331 pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 6 00332 pc_nDarray_rep->nDarray_rank = rank_of_nDarray; //rank_of_nDarray; 00333 00334 // in the case of nDarray_rank=0 add one to get right thing from the 00335 // operator new 00336 int one_or0 = 0; 00337 if(!pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00338 pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00339 // dimensions 00340 00341 pc_nDarray_rep->total_numb = 1; 00342 for( int idim = 0 ; idim < pc_nDarray_rep->nDarray_rank ; idim++ ) 00343 { 00344 pc_nDarray_rep->dim[idim] = pdim[idim]; 00345 pc_nDarray_rep->total_numb *= pc_nDarray_rep->dim[idim]; 00346 } 00347 00348 // allocate memory for the actual nDarray as nDarray 00349 pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00350 if (!pc_nDarray_rep->pd_nDdata) 00351 { 00352 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 00353 ::exit(1); 00354 } 00355 00356 pc_nDarray_rep->n = 1; // so far, there's one reference 00357 00358 switch(pc_nDarray_rep->nDarray_rank) 00359 { 00360 case 0: 00361 { 00362 ::printf("\a\n Unit nDarray of rank 0 ???\n"); 00363 break; 00364 } 00365 00366 case 1: 00367 { 00368 ::printf("\a\n Unit nDarray of rank 1 ???\n"); 00369 break; 00370 } 00371 00372 case 2: // Kronecker delta 00373 00374 if ( flag[0] == 'I' ) 00375 { 00376 for ( int i2=1 ; i2<=pc_nDarray_rep->dim[0] ; i2++ ) 00377 { 00378 for ( int j2=1 ; j2<=pc_nDarray_rep->dim[1] ; j2++ ) 00379 { 00380 val(i2,j2) = (i2 == j2 ? 1 : 0); 00381 } 00382 } 00383 break; 00384 00385 } 00386 00387 case 3: // Levi - Civita permutation BJtensor 00388 { 00389 for ( int i3=1 ; i3<=pc_nDarray_rep->dim[0] ; i3++ ) 00390 { 00391 for ( int j3=1 ; j3<=pc_nDarray_rep->dim[1] ; j3++ ) 00392 { 00393 for ( int k3=1 ; k3<=pc_nDarray_rep->dim[2] ; k3++ ) 00394 { 00395 if ( i3==1 && j3==2 && k3==3 || 00396 i3==2 && j3==3 && k3==1 || 00397 i3==3 && j3==1 && k3==2 ) 00398 { 00399 val(i3,j3,k3) = 1.0; 00400 } 00401 else if ( i3==3 && j3==2 && k3==1 || 00402 i3==2 && j3==1 && k3==3 || 00403 i3==1 && j3==3 && k3==2 ) 00404 { 00405 val(i3,j3,k3) = -1.0; 00406 } 00407 else val(i3,j3,k3) = 0.0; 00408 } 00409 } 00410 } 00411 break; 00412 } 00413 00414 // case 4: 00415 // for ( int i4=1 ; i4<=pc_nDarray_rep->dim[0] ; i4++ ) 00416 // for ( int j4=1 ; j4<=pc_nDarray_rep->dim[1] ; j4++ ) 00417 // for ( int k4=1 ; k4<=pc_nDarray_rep->dim[2] ; k4++ ) 00418 // for ( int l4=1 ; l4<=pc_nDarray_rep->dim[3] ; l4++ ) 00419 // val(i4,j4,k4,l4) = (((i4==k4) && (j4==l4)) ? 1 : 0); 00420 // break; 00421 } 00422 00423 00424 } 00425 00426 00427 //############################################################################## 00428 nDarray::nDarray(const nDarray & x) 00429 { 00430 x.pc_nDarray_rep->n++; // we're adding another reference. 00431 pc_nDarray_rep = x.pc_nDarray_rep; // point to the new nDarray_rep. 00432 } 00433 00434 00435 00436 //############################################################################## 00437 nDarray::~nDarray() 00438 { 00439 if (--pc_nDarray_rep->n == 0) // if reference count goes to 0 00440 { 00441 // DEallocate memory of the actual nDarray 00442 // delete [pc_nDarray_rep->pc_nDarray_rep->total_numb] pc_nDarray_rep->pd_nDdata; 00443 // see ELLIS & STROUSTRUP $18.3 00444 // and note on the p.65($5.3.4) 00445 // and the page 276 ($12.4) 00446 delete [] pc_nDarray_rep->pd_nDdata; 00447 delete [] pc_nDarray_rep->dim; 00448 delete pc_nDarray_rep; 00449 } 00450 } 00451 00452 00453 //############################################################################## 00454 // use "from" and initialize already allocated BJtensor from "from" values 00455 void nDarray::Initialize( const nDarray & from ) 00456 { 00457 // copy only data because everything else has already been defined 00458 // WATCH OUT IT HAS TO BE DEFINED BEFORE THIS FUNCTIONS IS CALLED 00459 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00460 this->pc_nDarray_rep->pd_nDdata[i] = from.pc_nDarray_rep->pd_nDdata[i] ; 00461 } 00462 //############################################################################## 00463 // // reset data to "value" 00464 void nDarray::Reset_to( double value ) // reset data to "value" 00465 { 00466 // set only data because everything else has already been defined 00467 // WATCH OUT IT HAS TO BE DEFINED BEFORE THIS FUNCTIONS IS CALLED 00468 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00469 this->pc_nDarray_rep->pd_nDdata[i] = value; 00470 } 00471 //############################################################################## 00472 // use "from" and initialize AND allocate BJtensor from "from" values 00473 // it should work on all the things that are missed by default constructor 00474 //(int rank_of_nDarray, const int *pdim, double *values) 00475 void nDarray::Initialize_all( const nDarray & from ) 00476 { 00477 // create the structure: 00478 // pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded //IN DEFAULT 00479 this->pc_nDarray_rep->nDarray_rank = from.rank(); //rank_of_nDarray; 00480 // in the case of nDarray_rank=0 add one to get right thing from the 00481 // operator new 00482 int one_or0 = 0; 00483 if(!this->pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00484 delete [] this->pc_nDarray_rep->dim; // get rid of old allocated memory 00485 this->pc_nDarray_rep->dim = new int[pc_nDarray_rep->nDarray_rank+one_or0];// array for 00486 // dimensions 00487 this->pc_nDarray_rep->total_numb = 1; 00488 for( int idim = 0 ; idim < this->pc_nDarray_rep->nDarray_rank ; idim++ ) 00489 { 00490 this->pc_nDarray_rep->dim[idim] = from.dim()[idim]; // fill dims from from!! 00491 this->pc_nDarray_rep->total_numb *= pc_nDarray_rep->dim[idim]; // find total number 00492 } 00493 delete [] this->pc_nDarray_rep->pd_nDdata; // get rid of old allocated memory 00494 // allocate memory for the actual nDarray as nDarray 00495 this->pc_nDarray_rep->pd_nDdata = new double [(size_t)pc_nDarray_rep->total_numb]; 00496 if (!this->pc_nDarray_rep->pd_nDdata) 00497 { 00498 ::fprintf(stderr,"\a\nInsufficient memory for array in Initialize_all \n"); 00499 ::exit(1); 00500 } 00501 this->pc_nDarray_rep->n = 1; // so far, there's one reference 00502 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 00503 this->pc_nDarray_rep->pd_nDdata[i] = from.pc_nDarray_rep->pd_nDdata[i]; 00504 00505 } 00506 00507 //############################################################################## 00508 nDarray nDarray::deep_copy() 00509 { 00510 // create the structure: 00511 nDarray temp; 00512 temp.pc_nDarray_rep->nDarray_rank = this->pc_nDarray_rep->nDarray_rank; 00513 00514 // in the case of nDarray_rank=0 add one to get right thing from the 00515 // operator new 00516 int one_or0 = 0; 00517 if(!temp.pc_nDarray_rep->nDarray_rank) one_or0 = 1; 00518 delete [] temp.pc_nDarray_rep->dim; //delete default value 00519 temp.pc_nDarray_rep->dim = new int[temp.pc_nDarray_rep->nDarray_rank+one_or0];// array for 00520 // dimensions 00521 00522 for( int idim = 0 ; idim < temp.pc_nDarray_rep->nDarray_rank ; idim++ ) 00523 { 00524 temp.pc_nDarray_rep->dim[idim] = this->pc_nDarray_rep->dim[idim] ; 00525 } 00526 temp.pc_nDarray_rep->total_numb = this->pc_nDarray_rep->total_numb; 00527 00528 delete [] temp.pc_nDarray_rep->pd_nDdata;//delete default value 00529 temp.pc_nDarray_rep->pd_nDdata = new double [temp.pc_nDarray_rep->total_numb]; 00530 if (!temp.pc_nDarray_rep->pd_nDdata) 00531 { 00532 ::fprintf(stderr,"\a\nInsufficient memory for array in deep_copy\n"); 00533 ::exit(1); 00534 } 00535 00536 temp.pc_nDarray_rep->n = 1; // so far, there's one reference 00537 00538 for ( int i=0 ; i<temp.pc_nDarray_rep->total_numb ; i++ ) 00539 temp.pc_nDarray_rep->pd_nDdata[i] = this->pc_nDarray_rep->pd_nDdata[i] ; 00540 00541 return temp; 00542 } 00543 00544 00545 //############################################################################## 00546 nDarray& nDarray::operator=(const nDarray & rval) 00547 { 00548 rval.pc_nDarray_rep->n++; // tell the rval it has another reference 00549 00550 // /* It is important to increment the reference_counter in the new 00551 // nDarray before decrementing the reference_counter in the 00552 // old nDarray_rep to ensure proper operation when assigning a 00553 // nDarray_rep to itself ( after ARKoenig JOOP May/June '90 ) */ 00554 00555 // clean up current value; 00556 if(--pc_nDarray_rep->n == 0) // if nobody else is referencing us. 00557 { 00558 // DEallocate memory of the actual nDarray 00559 // delete [pc_nDarray_rep->pc_nDarray_rep->total_numb] pc_nDarray_rep->pd_nDdata; 00560 // see ELLIS & STROUSTRUP $18.3 00561 // and note on the p.65($5.3.4) 00562 delete [] pc_nDarray_rep->pd_nDdata; 00563 delete [] pc_nDarray_rep->dim; 00564 delete pc_nDarray_rep; 00565 } 00566 00567 // connect to new value 00568 pc_nDarray_rep = rval.pc_nDarray_rep; // point at the rval nDarray_rep 00569 return *this; 00570 } 00571 00572 00573 //############################################################################## 00574 double & nDarray::val(int subscript, ...) 00575 { 00576 // if scalar get back 00577 if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00578 // for all others procede 00579 va_list p_arg; 00580 va_start(p_arg, subscript); // initialize p_arg 00581 00582 long int first = 0; 00583 long int second = 0; 00584 00585 first = subscript; // first index 00586 long int where = first - 1; 00587 for ( int Dcount=1 ; Dcount<=pc_nDarray_rep->nDarray_rank-1 ; Dcount++ ) 00588 { // for all dimensions less then 1 this will be skipped 00589 second = va_arg(p_arg, int); // next 00590 where = where*pc_nDarray_rep->dim[Dcount]+second - 1; 00591 first = second; 00592 } 00593 va_end(p_arg); 00594 //::printf("*w=%2ld ",where); 00595 double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00596 return (*p_value); 00597 } 00598 // ..JB..//############################################################################## 00599 // ..JB..// another overloading of operator() . . . // overloaded for ONE argument 00600 // ..JB..double & nDarray::val(int first) 00601 // ..JB.. { 00602 // ..JB.. long int where = first - 1; 00603 // ..JB.. 00604 // ..JB..//::printf(" w=%ld ",where); 00605 // ..JB.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00606 // ..JB.. return (*p_value); 00607 // ..JB.. } 00608 // ..JB.. 00609 // ..JB..//############################################################################## 00610 // ..JB..// another overloading of operator() . . . // overloaded for TWO arguments 00611 // ..JB..double & nDarray::val(int first, int second) 00612 // ..JB.. { 00613 // ..JB.. long int where = first - 1; 00614 // ..JB.. where = where*pc_nDarray_rep->dim[1]+second - 1; 00615 // ..JB.. 00616 // ..JB..//::printf(" w=%ld ",where); 00617 // ..JB.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00618 // ..JB.. return (*p_value); 00619 // ..JB.. } 00620 // ..JB.. 00621 // ..JB..//############################################################################## 00622 // ..JB..// another overloading of operator() . . . // overloaded for THREE arguments 00623 // ..JB..double & nDarray::val(int first, int second, int third ) 00624 // ..JB.. { 00625 // ..JB.. long int where = first - 1; 00626 // ..JB.. where = where*pc_nDarray_rep->dim[1]+second - 1; 00627 // ..JB.. where = where*pc_nDarray_rep->dim[2]+third - 1; 00628 // ..JB.. 00629 // ..JB..//::printf(" w=%ld ",where); 00630 // ..JB.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00631 // ..JB.. return (*p_value); 00632 // ..JB.. } 00633 // ..JB.. 00634 // ..JB..//############################################################################## 00635 // ..JB..// another overloading of operator() . . . // overloaded for FOUR arguments 00636 // ..JB..double & nDarray::val(int first, int second, 00637 // ..JB.. int third, int fourth) 00638 // ..JB.. { 00639 // ..JB.. if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00640 // ..JB.. 00641 // ..JB.. long int where = first - 1; 00642 // ..JB.. 00643 // ..JB.. if(pc_nDarray_rep->nDarray_rank==2) 00644 // ..JB.. { 00645 // ..JB.. where = where*pc_nDarray_rep->dim[1]+second - 1; 00646 // ..JB.. } 00647 // ..JB.. 00648 // ..JB.. if(pc_nDarray_rep->nDarray_rank==3) 00649 // ..JB.. { 00650 // ..JB.. where = where*pc_nDarray_rep->dim[2]+third - 1; 00651 // ..JB.. } 00652 // ..JB.. 00653 // ..JB.. if(pc_nDarray_rep->nDarray_rank==4) 00654 // ..JB.. { 00655 // ..JB.. where = where*pc_nDarray_rep->dim[3]+fourth - 1; 00656 // ..JB.. } 00657 // ..JB.. 00658 // ..JB..//::printf(" w=%ld ",where); 00659 // ..JB.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00660 // ..JB.. return (*p_value); 00661 // ..JB.. } 00662 // ..JB.. 00663 // ..JB.. 00664 00665 00666 00667 //############################################################################## 00668 // another overloading of operator() . . . // overloaded for FOUR arguments 00669 double & nDarray::val4(int first, int second, 00670 int third, int fourth) 00671 { 00672 //JB//JB long int where = (((first-1)*pc_nDarray_rep->dim[1]+second-1)*pc_nDarray_rep->dim[2]+third-1)*pc_nDarray_rep->dim[3]+fourth-1; 00673 //JB double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00674 //JB// ..JB.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00675 //JB return (*p_value); 00676 //JB// return (*pc_nDarray_rep->pd_nDdata + (size_t)(((first-1)*pc_nDarray_rep->dim[1]+second-1)*pc_nDarray_rep->dim[2]+third-1)*pc_nDarray_rep->dim[3]+fourth-1); 00677 00678 00679 if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00680 00681 long int where = first - 1; 00682 00683 if(pc_nDarray_rep->nDarray_rank==2) 00684 { 00685 where = where*pc_nDarray_rep->dim[1]+second - 1; 00686 } 00687 00688 if(pc_nDarray_rep->nDarray_rank==3) 00689 { 00690 where = where*pc_nDarray_rep->dim[2]+third - 1; 00691 } 00692 00693 if(pc_nDarray_rep->nDarray_rank==4) 00694 { 00695 where = where*pc_nDarray_rep->dim[3]+fourth - 1; 00696 } 00697 00698 //::printf(" w=%ld ",where); 00699 double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00700 return (*p_value); 00701 00702 00703 00704 } 00705 00706 00707 00708 00709 //..//############################################################################## 00710 //..// another overloading of operator() . . . 00711 //..// overloaded for more than FOUR arguments 00712 //..double & nDarray::val(int first, int second, 00713 //.. int third, int fourth, 00714 //.. int subscript, ...) 00715 //.. { 00716 //..// ako je skalar tojest nulti red onda o'ma nazad 00717 //.. if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00718 //..// za bilo koji veci red ajd' dalje 00719 //.. va_list p_arg; 00720 //.. va_start(p_arg, subscript); // initialize p_arg 00721 //.. 00722 //.. int prvi = subscript; 00723 //.. int drugi = 0; 00724 //.. 00725 //.. long int where = first - 1; 00726 //.. where = where*pc_nDarray_rep->dim[1]+second - 1; 00727 //.. where = where*pc_nDarray_rep->dim[2]+third - 1; 00728 //.. where = where*pc_nDarray_rep->dim[3]+fourth - 1; 00729 //.. 00730 //..// dalje iza cetvrtog 00731 //.. where = where*pc_nDarray_rep->dim[4]+prvi - 1; 00732 //.. 00733 //.. 00734 //.. for ( int Dcount=5 ; Dcount<=pc_nDarray_rep->nDarray_rank-1 ; Dcount++ ) 00735 //.. { // ovo ce ustvari biti preskoceno za sve dimenzije manje od 2 00736 //.. drugi = va_arg(p_arg, int); // sledeci 00737 //.. where = where*pc_nDarray_rep->dim[Dcount]+drugi - 1; 00738 //.. prvi = drugi; 00739 //.. } 00740 //.. va_end(p_arg); 00741 //.. 00742 //..//::printf(" w=%ld ",where); 00743 //.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00744 //.. return (*p_value); 00745 //.. } 00746 00747 00748 00749 //############################################################################## 00750 double nDarray::cval(int subscript, ...) const 00751 { 00752 // if scalar get back 00753 if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00754 // for all others procede 00755 va_list p_arg; 00756 va_start(p_arg, subscript); // initialize p_arg 00757 00758 int first = 0; 00759 int second = 0; 00760 00761 first = subscript; // first indeks 00762 long int where = first - 1; 00763 for ( int Dcount=1 ; Dcount<=pc_nDarray_rep->nDarray_rank-1 ; Dcount++ ) 00764 { // for all dimensions less then 2 this will be skipped 00765 second = va_arg(p_arg, int); // next 00766 where = where*pc_nDarray_rep->dim[Dcount]+second - 1; 00767 first = second; 00768 } 00769 va_end(p_arg); 00770 00771 //::printf("*w=%2ld ",where); 00772 double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00773 return (*p_value); 00774 } 00775 00776 00777 //..//############################################################################## 00778 //..// another overloading of operator() . . . 00779 //..double nDarray::operator()(int subscript, ...) const 00780 //.. { 00781 //..// ako je skalar tojest nulti red onda o'ma nazad 00782 //.. if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00783 //..// za bilo koji veci red ajd' dalje 00784 //.. va_list p_arg; 00785 //.. va_start(p_arg, subscript); // initialize p_arg 00786 //.. 00787 //.. int prvi = 0; 00788 //.. int drugi = 0; 00789 //.. 00790 //.. prvi = subscript; // prvi indeks 00791 //.. long int where = prvi - 1; 00792 //.. for ( int Dcount=1 ; Dcount<=pc_nDarray_rep->nDarray_rank-1 ; Dcount++ ) 00793 //.. { // ovo ce ustvari biti preskoceno za sve dimenzije manje od 2 00794 //.. drugi = va_arg(p_arg, int); // sledeci 00795 //.. where = where*pc_nDarray_rep->dim[Dcount]+drugi - 1; 00796 //.. prvi = drugi; 00797 //.. } 00798 //.. va_end(p_arg); 00799 //.. 00800 //..//::printf(" w=%ld ",where); 00801 //.. double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00802 //.. return (*p_value); 00803 //.. } 00804 //.. 00805 //.. 00806 //.. 00807 00808 //@@@@@ 00809 //@@@@@ //############################################################################## 00810 //@@@@@ // another overloading of operator() . . . // overloaded for ONE argument 00811 //@@@@@ double nDarray::operator()(int first) const 00812 //@@@@@ { 00813 //@@@@@ long int where = first - 1; 00814 //@@@@@ 00815 //@@@@@ //::printf(" w=%ld ",where); 00816 //@@@@@ double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00817 //@@@@@ return (*p_value); 00818 //@@@@@ } 00819 //@@@@@ 00820 //@@@@@ //############################################################################## 00821 //@@@@@ // another overloading of operator() . . . // overloaded for TWO arguments 00822 //@@@@@ double nDarray::operator()(int first, int second) const 00823 //@@@@@ { 00824 //@@@@@ long int where = first - 1; 00825 //@@@@@ where = where*pc_nDarray_rep->dim[1]+second - 1; 00826 //@@@@@ 00827 //@@@@@ //::printf(" w=%ld ",where); 00828 //@@@@@ double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00829 //@@@@@ return (*p_value); 00830 //@@@@@ } 00831 //@@@@@ 00832 //@@@@@ //############################################################################## 00833 //@@@@@ // another overloading of operator() . . . // overloaded for THREE arguments 00834 //@@@@@ double nDarray::operator()(int first, int second, int third ) const 00835 //@@@@@ { 00836 //@@@@@ long int where = first - 1; 00837 //@@@@@ where = where*pc_nDarray_rep->dim[1]+second - 1; 00838 //@@@@@ where = where*pc_nDarray_rep->dim[2]+third - 1; 00839 //@@@@@ 00840 //@@@@@ //::printf(" w=%ld ",where); 00841 //@@@@@ double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00842 //@@@@@ return (*p_value); 00843 //@@@@@ } 00844 //@@@@@ 00845 //@@@@@ //############################################################################## 00846 //@@@@@ // another overloading of operator() . . . // overloaded for FOUR arguments 00847 //@@@@@ double nDarray::operator()(int first, int second, 00848 //@@@@@ int third, int fourth) const 00849 //@@@@@ { 00850 //@@@@@ long int where = first - 1; 00851 //@@@@@ where = where*pc_nDarray_rep->dim[1]+second - 1; 00852 //@@@@@ where = where*pc_nDarray_rep->dim[2]+third - 1; 00853 //@@@@@ where = where*pc_nDarray_rep->dim[3]+fourth - 1; 00854 //@@@@@ 00855 //@@@@@ //::printf(" w=%ld ",where); 00856 //@@@@@ double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00857 //@@@@@ return (*p_value); 00858 //@@@@@ } 00859 //@@@@@ 00860 //@@@@@ 00861 //@@@@@ 00862 //@@@@@ //############################################################################## 00863 //@@@@@ // another overloading of operator() . . . 00864 //@@@@@ // overloaded for more than FOUR arguments 00865 //@@@@@ double nDarray::operator()(int first, int second, 00866 //@@@@@ int third, int fourth, 00867 //@@@@@ int subscript, ...) const 00868 //@@@@@ { 00869 //@@@@@ // ako je skalar tojest nulti red onda o'ma nazad 00870 //@@@@@ if(pc_nDarray_rep->nDarray_rank==0) return (*pc_nDarray_rep->pd_nDdata); 00871 //@@@@@ // za bilo koji veci red ajd' dalje 00872 //@@@@@ va_list p_arg; 00873 //@@@@@ va_start(p_arg, subscript); // initialize p_arg 00874 //@@@@@ 00875 //@@@@@ int prvi = subscript; 00876 //@@@@@ int drugi = 0; 00877 //@@@@@ 00878 //@@@@@ long int where = first - 1; 00879 //@@@@@ where = where*pc_nDarray_rep->dim[1]+second - 1; 00880 //@@@@@ where = where*pc_nDarray_rep->dim[2]+third - 1; 00881 //@@@@@ where = where*pc_nDarray_rep->dim[3]+fourth - 1; 00882 //@@@@@ 00883 //@@@@@ // dalje iza cetvrtog 00884 //@@@@@ where = where*pc_nDarray_rep->dim[4]+prvi - 1; 00885 //@@@@@ 00886 //@@@@@ 00887 //@@@@@ for ( int Dcount=5 ; Dcount<=pc_nDarray_rep->nDarray_rank-1 ; Dcount++ ) 00888 //@@@@@ { // ovo ce ustvari biti preskoceno za sve dimenzije manje od 2 00889 //@@@@@ drugi = va_arg(p_arg, int); // sledeci 00890 //@@@@@ where = where*pc_nDarray_rep->dim[Dcount]+drugi - 1; 00891 //@@@@@ prvi = drugi; 00892 //@@@@@ } 00893 //@@@@@ va_end(p_arg); 00894 //@@@@@ 00895 //@@@@@ //::printf(" w=%ld ",where); 00896 //@@@@@ double *p_value = pc_nDarray_rep->pd_nDdata + (size_t)where; 00897 //@@@@@ return (*p_value); 00898 //@@@@@ } 00899 //@@@@@ 00900 //@@@@@ 00901 //@@@@@ 00902 00903 00904 00905 //++//############################################################################## 00906 //++// nDarray addition 00907 //++nDarray nDarray::operator+( nDarray & rval) 00908 //++ { 00909 //++ int this_rank_of_nDarray = this->pc_nDarray_rep->nDarray_rank; 00910 //++ int rval_rank_of_nDarray = rval.pc_nDarray_rep->nDarray_rank; 00911 //++ 00912 //++ if(this_rank_of_nDarray != rval_rank_of_nDarray) 00913 //++ { 00914 //++ ::printf("\a\nnDarrays of different ranks: addition not possible\n"); 00915 //++ ::exit ( 1 ); 00916 //++ } 00917 //++ 00918 //++ for ( int i=0 ; i<this_rank_of_nDarray ; i++ ) 00919 //++ if (this->pc_nDarray_rep->dim[i] != rval.pc_nDarray_rep->dim[i] ) 00920 //++ { 00921 //++ ::fprintf(stderr,"\a\nDimension discrepancy in operator+\n", 00922 //++ "this->pc_nDarray_rep->dim[%d]=%d\n", 00923 //++ "arg.pc_nDarray_rep->dim[%d]=%d\n", 00924 //++ i,this->pc_nDarray_rep->dim[i], 00925 //++ i,rval.pc_nDarray_rep->dim[i]); 00926 //++ ::exit(1); 00927 //++ } 00928 //++// construct nDarray using the same control numbers as for the 00929 //++// original one . 00930 //++ nDarray add(pc_nDarray_rep->nDarray_rank, pc_nDarray_rep->dim, 0.0); 00931 //++ 00932 //++ switch(pc_nDarray_rep->nDarray_rank) 00933 //++ { 00934 //++ case 0: 00935 //++ { 00936 //++ add.val(1) = val(1) + rval.val(1); 00937 //++ break; 00938 //++ } 00939 //++ 00940 //++ case 1: 00941 //++ { 00942 //++ for ( int i1=1 ; i1<=this->pc_nDarray_rep->dim[0] ; i1++ ) 00943 //++ { 00944 //++ add.val(i1) = val(i1) + rval.val(i1); 00945 //++ } 00946 //++ break; 00947 //++ } 00948 //++ 00949 //++ case 2: 00950 //++ { 00951 //++ for ( int i2=1 ; i2<=this->pc_nDarray_rep->dim[0] ; i2++ ) 00952 //++ { 00953 //++ for ( int j2=1 ; j2<=this->pc_nDarray_rep->dim[1] ; j2++ ) 00954 //++ { 00955 //++ add.val(i2, j2) = val(i2, j2) + rval.val(i2, j2); 00956 //++ } 00957 //++ } 00958 //++ break; 00959 //++ } 00960 //++ 00961 //++ case 3: 00962 //++ { 00963 //++ for ( int i3=1 ; i3<=this->pc_nDarray_rep->dim[0] ; i3++ ) 00964 //++ { 00965 //++ for ( int j3=1 ; j3<=this->pc_nDarray_rep->dim[1] ; j3++ ) 00966 //++ { 00967 //++ for ( int k3=1 ; k3<=this->pc_nDarray_rep->dim[2] ; k3++ ) 00968 //++ { 00969 //++ add.val(i3, j3, k3) = val(i3, j3, k3) + rval.val(i3, j3, k3); 00970 //++ } 00971 //++ } 00972 //++ } 00973 //++ break; 00974 //++ } 00975 //++ 00976 //++ case 4: 00977 //++ { 00978 //++ for ( int i4=1 ; i4<=this->pc_nDarray_rep->dim[0] ; i4++ ) 00979 //++ { 00980 //++ for ( int j4=1 ; j4<=this->pc_nDarray_rep->dim[1] ; j4++ ) 00981 //++ { 00982 //++ for ( int k4=1 ; k4<=this->pc_nDarray_rep->dim[2] ; k4++ ) 00983 //++ { 00984 //++ for ( int l4=1 ; l4<=this->pc_nDarray_rep->dim[3] ; l4++ ) 00985 //++ { 00986 //++ add.val(i4,j4,k4,l4)=val(i4,j4,k4,l4)+rval.val(i4,j4,k4,l4); 00987 //++ } 00988 //++ } 00989 //++ } 00990 //++ } 00991 //++ break; 00992 //++ } 00993 //++ } 00994 //++ 00995 //++ return add; 00996 //++ } 00997 00998 //############################################################################## 00999 // nDarray addition 01000 nDarray& nDarray::operator+=(const nDarray & rval) 01001 { 01002 int this_rank_of_nDarray = this->pc_nDarray_rep->nDarray_rank; 01003 int rval_rank_of_nDarray = rval.pc_nDarray_rep->nDarray_rank; 01004 01005 if(this_rank_of_nDarray != rval_rank_of_nDarray) 01006 { 01007 ::printf("\a\nnDarrays of different ranks: += not possible\n"); 01008 ::exit ( 1 ); 01009 } 01010 01011 int i = 0; 01012 for ( i=0 ; i<this_rank_of_nDarray ; i++ ) 01013 if (this->pc_nDarray_rep->dim[i] != rval.pc_nDarray_rep->dim[i] ) 01014 { 01015 ::fprintf(stderr,"\a\nDimension discrepancy in operator+=\n\ 01016 this->pc_nDarray_rep->dim[%d]=%d\n\ 01017 arg.pc_nDarray_rep->dim[%d]=%d\n", 01018 i,this->pc_nDarray_rep->dim[i], 01019 i,rval.pc_nDarray_rep->dim[i]); 01020 ::exit(1); 01021 } 01022 // Copy *this if necessary 01023 if ( this->pc_nDarray_rep->n > 1 )// see ARK in JOOP may/june '90 01024 { // "Letter From a Newcomer" 01025 //.............................................................................. 01026 // create the structure: 01027 nDarray_rep * New_pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 01028 New_pc_nDarray_rep->nDarray_rank = this->pc_nDarray_rep->nDarray_rank; 01029 // in the case of nDarray_rank=0 add one to get right thing from the 01030 // operator new 01031 int one_or0 = 0; 01032 if(!New_pc_nDarray_rep->nDarray_rank) one_or0 = 1; 01033 New_pc_nDarray_rep->dim = new int[New_pc_nDarray_rep->nDarray_rank+one_or0]; 01034 // array for dimensions 01035 New_pc_nDarray_rep->total_numb = 1; 01036 for( int idim = 0 ; idim < New_pc_nDarray_rep->nDarray_rank ; idim++ ) 01037 { 01038 New_pc_nDarray_rep->dim[idim] = this->pc_nDarray_rep->dim[idim]; 01039 New_pc_nDarray_rep->total_numb *= New_pc_nDarray_rep->dim[idim]; 01040 } 01041 // allocate memory for the actual nDarray as nDarray 01042 New_pc_nDarray_rep->pd_nDdata = new double [(size_t)New_pc_nDarray_rep->total_numb]; 01043 if (!New_pc_nDarray_rep->pd_nDdata) 01044 { 01045 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 01046 ::exit(1); 01047 } 01048 New_pc_nDarray_rep->n = 1; // so far, there's one reference 01049 for ( i=0 ; i<New_pc_nDarray_rep->total_numb ; i++ ) 01050 New_pc_nDarray_rep->pd_nDdata[i] = this->pc_nDarray_rep->pd_nDdata[i]; 01051 //......... 01052 this->pc_nDarray_rep->n--; 01053 this->pc_nDarray_rep = New_pc_nDarray_rep; 01054 //.............................................................................. 01055 } 01056 // it appears that I can add this two nDarrays just as a simple BJvectors: 01057 for (int j=0 ; j<this->pc_nDarray_rep->total_numb ; j++) 01058 this->pc_nDarray_rep->pd_nDdata[j] += rval.pc_nDarray_rep->pd_nDdata[j]; 01059 01060 return *this; 01061 } 01062 01063 01064 //############################################################################## 01065 // nDarray addition 01066 #ifndef _VC6 01067 nDarray operator+(const nDarray & lval, const nDarray & rval) 01068 { 01069 nDarray result(lval); 01070 result += rval; 01071 return result; 01072 } 01073 #else 01074 nDarray nDarray::operator+(const nDarray & rval) const 01075 { 01076 nDarray result(*this); 01077 result += rval; 01078 return result; 01079 } 01080 #endif 01081 01082 01083 01084 //############################################################################## 01085 // scalar addition 01086 nDarray nDarray::operator+( double rval) 01087 { 01088 // construct nDarray using the same control numbers as for the 01089 // original one. 01090 nDarray add(pc_nDarray_rep->nDarray_rank, pc_nDarray_rep->dim, 0.0); 01091 switch(pc_nDarray_rep->nDarray_rank) 01092 { 01093 case 0: 01094 { 01095 add.val(1) = val(1) + rval; 01096 break; 01097 } 01098 case 1: 01099 { 01100 for ( int i1=1 ; i1<=this->pc_nDarray_rep->dim[0] ; i1++ ) 01101 { 01102 add.val(i1) = val(i1) + rval; 01103 } 01104 break; 01105 } 01106 case 2: 01107 { 01108 for ( int i2=1 ; i2<=this->pc_nDarray_rep->dim[0] ; i2++ ) 01109 { 01110 for ( int j2=1 ; j2<=this->pc_nDarray_rep->dim[1] ; j2++ ) 01111 { 01112 add.val(i2, j2) = val(i2, j2) + rval; 01113 } 01114 } 01115 break; 01116 } 01117 case 3: 01118 { 01119 for ( int i3=1 ; i3<=this->pc_nDarray_rep->dim[0] ; i3++ ) 01120 { 01121 for ( int j3=1 ; j3<=this->pc_nDarray_rep->dim[1] ; j3++ ) 01122 { 01123 for ( int k3=1 ; k3<=this->pc_nDarray_rep->dim[2] ; k3++ ) 01124 { 01125 add.val(i3, j3, k3) = val(i3, j3, k3) + rval; 01126 } 01127 } 01128 } 01129 break; 01130 } 01131 case 4: 01132 { 01133 for ( int i4=1 ; i4<=this->pc_nDarray_rep->dim[0] ; i4++ ) 01134 { 01135 for ( int j4=1 ; j4<=this->pc_nDarray_rep->dim[1] ; j4++ ) 01136 { 01137 for ( int k4=1 ; k4<=this->pc_nDarray_rep->dim[2] ; k4++ ) 01138 { 01139 for ( int l4=1 ; l4<=this->pc_nDarray_rep->dim[3] ; l4++ ) 01140 { 01141 add.val(i4,j4,k4,l4)=val(i4,j4,k4,l4)+rval; 01142 } 01143 } 01144 } 01145 } 01146 break; 01147 } 01148 } 01149 return add; 01150 } 01151 01152 //############################################################################## 01153 // scalar multiplication 01154 nDarray nDarray::operator*( const double rval) const 01155 { 01156 // construct nDarray using the same control numbers as for the 01157 // original one. 01158 nDarray mult(pc_nDarray_rep->nDarray_rank, pc_nDarray_rep->dim, 0.0); 01159 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 01160 mult.pc_nDarray_rep->pd_nDdata[i] = this->pc_nDarray_rep->pd_nDdata[i]*rval; 01161 01162 return mult; 01163 } 01164 01165 // nDarray operator*( double lval, nDarray & rval); // REVIEWER global 01166 //############################################################################## 01167 // scalar multiplication 01168 nDarray operator*( const double lval, const nDarray & rval) 01169 { 01170 return rval*lval; 01171 } 01172 //############################################################################## 01173 // nDarray substraction 01174 nDarray& nDarray::operator-=(const nDarray & rval) 01175 { 01176 int this_rank_of_nDarray = this->pc_nDarray_rep->nDarray_rank; 01177 int rval_rank_of_nDarray = rval.pc_nDarray_rep->nDarray_rank; 01178 if(this_rank_of_nDarray != rval_rank_of_nDarray) 01179 { 01180 ::printf("\a\nnDarrays of different ranks: -= not possible\n"); 01181 ::exit ( 1 ); 01182 } 01183 01184 int i = 0; 01185 for ( i=0 ; i<this_rank_of_nDarray ; i++ ) 01186 if (this->pc_nDarray_rep->dim[i] != rval.pc_nDarray_rep->dim[i] ) 01187 { 01188 ::fprintf(stderr,"\a\nDimension discrepancy in operator+\n\ 01189 this->pc_nDarray_rep->dim[%d]=%d\n\ 01190 arg.pc_nDarray_rep->dim[%d]=%d\n", 01191 i,this->pc_nDarray_rep->dim[i], 01192 i,rval.pc_nDarray_rep->dim[i]); 01193 ::exit(1); 01194 } 01195 // Copy *this if necessary 01196 if ( this->pc_nDarray_rep->n > 1 )// see ARK in JOOP may/june '90 01197 { // "Letter From a Newcomer" 01198 //.............................................................................. 01199 // create the structure: 01200 nDarray_rep * New_pc_nDarray_rep = new nDarray_rep; // this 'new' is overloaded 01201 New_pc_nDarray_rep->nDarray_rank = this->pc_nDarray_rep->nDarray_rank; 01202 // in the case of nDarray_rank=0 add one to get right thing from the 01203 // operator new 01204 int one_or0 = 0; 01205 if(!New_pc_nDarray_rep->nDarray_rank) one_or0 = 1; 01206 New_pc_nDarray_rep->dim = new int[New_pc_nDarray_rep->nDarray_rank+one_or0]; 01207 // array for dimensions 01208 New_pc_nDarray_rep->total_numb = 1; 01209 for( int idim = 0 ; idim < New_pc_nDarray_rep->nDarray_rank ; idim++ ) 01210 { 01211 New_pc_nDarray_rep->dim[idim] = this->pc_nDarray_rep->dim[idim]; 01212 New_pc_nDarray_rep->total_numb *= New_pc_nDarray_rep->dim[idim]; 01213 } 01214 // allocate memory for the actual nDarray as nDarray 01215 New_pc_nDarray_rep->pd_nDdata = new double [(size_t)New_pc_nDarray_rep->total_numb]; 01216 if (!New_pc_nDarray_rep->pd_nDdata) 01217 { 01218 ::fprintf(stderr,"\a\nInsufficient memory for array\n"); 01219 ::exit(1); 01220 } 01221 New_pc_nDarray_rep->n = 1; // so far, there's one reference 01222 for ( i=0 ; i<New_pc_nDarray_rep->total_numb ; i++ ) 01223 New_pc_nDarray_rep->pd_nDdata[i] = this->pc_nDarray_rep->pd_nDdata[i]; 01224 //......... 01225 this->pc_nDarray_rep->n--; 01226 this->pc_nDarray_rep = New_pc_nDarray_rep; 01227 //.............................................................................. 01228 } 01229 // it appears that I can add this two nDarrays just as a simple BJvectors: 01230 for (int j=0 ; j<this->pc_nDarray_rep->total_numb ; j++) 01231 this->pc_nDarray_rep->pd_nDdata[j] -= rval.pc_nDarray_rep->pd_nDdata[j]; 01232 return *this; 01233 } 01234 01235 01236 //############################################################################## 01237 // nDarray substraction 01238 nDarray operator-(const nDarray & lval, const nDarray & rval) 01239 { 01240 nDarray result(lval); 01241 result -= rval; 01242 return result; 01243 } 01244 01245 //++//############################################################################## 01246 //++// nDarray substraction 01247 //++nDarray nDarray::operator-( nDarray & rval) 01248 //++ { 01249 //++ int this_rank_of_nDarray = this->pc_nDarray_rep->nDarray_rank; 01250 //++ int rval_rank_of_nDarray = rval.pc_nDarray_rep->nDarray_rank; 01251 //++ 01252 //++ if(this_rank_of_nDarray != rval_rank_of_nDarray) 01253 //++ { 01254 //++ ::printf("\a\nnDarrays of different ranks:", 01255 //++ " substraction not possible\n"); 01256 //++ ::exit ( 1 ); 01257 //++ } 01258 //++ 01259 //++ for ( int i=0 ; i<this_rank_of_nDarray ; i++ ) 01260 //++ if (this->pc_nDarray_rep->dim[i] != rval.pc_nDarray_rep->dim[i] ) 01261 //++ { 01262 //++ ::fprintf(stderr,"\a\nDimension discrepancy in operator - \n", 01263 //++ "this->pc_nDarray_rep->dim[%d]=%d\n", 01264 //++ "arg.pc_nDarray_rep->dim[%d]=%d\n", 01265 //++ i,this->pc_nDarray_rep->dim[i], 01266 //++ i,rval.pc_nDarray_rep->dim[i]); 01267 //++ ::exit(1); 01268 //++ } 01269 //++// construct nDarray using the same control numbers as for the 01270 //++// original one. 01271 //++ nDarray sub(pc_nDarray_rep->nDarray_rank, pc_nDarray_rep->dim, 0.0); 01272 //++ 01273 //++ switch(pc_nDarray_rep->nDarray_rank) 01274 //++ { 01275 //++ case 0: 01276 //++ { 01277 //++ sub.val(1) = val(1) - rval.val(1); 01278 //++ break; 01279 //++ } 01280 //++ 01281 //++ case 1: 01282 //++ { 01283 //++ for ( int i1=1 ; i1<=this->pc_nDarray_rep->dim[0] ; i1++ ) 01284 //++ { 01285 //++ sub.val(i1) = val(i1) - rval.val(i1); 01286 //++ } 01287 //++ break; 01288 //++ } 01289 //++ 01290 //++ case 2: 01291 //++ { 01292 //++ for ( int i2=1 ; i2<=this->pc_nDarray_rep->dim[0] ; i2++ ) 01293 //++ { 01294 //++ for ( int j2=1 ; j2<=this->pc_nDarray_rep->dim[1] ; j2++ ) 01295 //++ { 01296 //++ sub.val(i2, j2) = val(i2, j2) - rval.val(i2, j2); 01297 //++ } 01298 //++ } 01299 //++ break; 01300 //++ } 01301 //++ 01302 //++ case 3: 01303 //++ { 01304 //++ for ( int i3=1 ; i3<=this->pc_nDarray_rep->dim[0] ; i3++ ) 01305 //++ { 01306 //++ for ( int j3=1 ; j3<=this->pc_nDarray_rep->dim[1] ; j3++ ) 01307 //++ { 01308 //++ for ( int k3=1 ; k3<=this->pc_nDarray_rep->dim[2] ; k3++ ) 01309 //++ { 01310 //++ sub.val(i3, j3, k3) = val(i3, j3, k3) - rval.val(i3, j3, k3); 01311 //++ } 01312 //++ } 01313 //++ } 01314 //++ break; 01315 //++ } 01316 //++ 01317 //++ case 4: 01318 //++ { 01319 //++ for ( int i4=1 ; i4<=this->pc_nDarray_rep->dim[0] ; i4++ ) 01320 //++ { 01321 //++ for ( int j4=1 ; j4<=this->pc_nDarray_rep->dim[1] ; j4++ ) 01322 //++ { 01323 //++ for ( int k4=1 ; k4<=this->pc_nDarray_rep->dim[2] ; k4++ ) 01324 //++ { 01325 //++ 01326 //++ 01327 //++ 01328 //++ 01329 //++ for ( int l4=1 ; l4<=this->pc_nDarray_rep->dim[3] ; l4++ ) 01330 //++ { 01331 //++ sub.val(i4,j4,k4,l4)=val(i4,j4,k4,l4)-rval.val(i4,j4,k4,l4); 01332 //++ } 01333 //++ } 01334 //++ } 01335 //++ } 01336 //++ break; 01337 //++ } 01338 //++ } 01339 //++ 01340 //++ return sub; 01341 //++ } 01342 //++ 01343 01344 //############################################################################## 01345 // scalar substraction 01346 nDarray nDarray::operator-( double rval) 01347 { 01348 // construct nDarray using the same control numbers as for the 01349 // original one. 01350 nDarray sub(pc_nDarray_rep->nDarray_rank, pc_nDarray_rep->dim, 0.0); 01351 01352 switch(pc_nDarray_rep->nDarray_rank) 01353 { 01354 case 0: 01355 { 01356 sub.val(1) = val(1) - rval; 01357 break; 01358 } 01359 01360 case 1: 01361 { 01362 for ( int i1=1 ; i1<=this->pc_nDarray_rep->dim[0] ; i1++ ) 01363 sub.val(i1) = val(i1) - rval; 01364 break; 01365 } 01366 01367 case 2: 01368 { 01369 for ( int i2=1 ; i2<=this->pc_nDarray_rep->dim[0] ; i2++ ) 01370 for ( int j2=1 ; j2<=this->pc_nDarray_rep->dim[1] ; j2++ ) 01371 sub.val(i2, j2) = val(i2, j2) - rval; 01372 break; 01373 } 01374 01375 case 3: 01376 { 01377 for ( int i3=1 ; i3<=this->pc_nDarray_rep->dim[0] ; i3++ ) 01378 for ( int j3=1 ; j3<=this->pc_nDarray_rep->dim[1] ; j3++ ) 01379 for ( int k3=1 ; k3<=this->pc_nDarray_rep->dim[2] ; k3++ ) 01380 sub.val(i3, j3, k3) = val(i3, j3, k3) - rval; 01381 break; 01382 } 01383 01384 case 4: 01385 { 01386 for ( int i4=1 ; i4<=this->pc_nDarray_rep->dim[0] ; i4++ ) 01387 for ( int j4=1 ; j4<=this->pc_nDarray_rep->dim[1] ; j4++ ) 01388 for ( int k4=1 ; k4<=this->pc_nDarray_rep->dim[2] ; k4++ ) 01389 for ( int l4=1 ; l4<=this->pc_nDarray_rep->dim[3] ; l4++ ) 01390 sub.val(i4,j4,k4,l4)=val(i4,j4,k4,l4)-rval; 01391 break; 01392 } 01393 } 01394 01395 return sub; 01396 } 01397 //############################################################################## 01398 //############################################################################## 01399 // unary minus 01400 nDarray nDarray::operator-( void ) 01401 { 01402 nDarray ret = *this *-1.0; 01403 01404 return ret; 01405 } 01406 01407 01408 //############################################################################# 01409 double nDarray::sum() const 01410 { 01411 double sum = 0; 01412 for ( int memb=0 ; memb<pc_nDarray_rep->total_numb ; memb++ ) 01413 sum += pc_nDarray_rep->pd_nDdata[memb]; 01414 return sum; 01415 } 01416 01417 //############################################################################## 01418 // trace function: trace of a second rank BJtensor 01419 // what about fourth ( 4th ) rank BJtensor trace or any other rank ?? 01420 double nDarray::trace() const 01421 { 01422 double tr = 0.0; 01423 // ovaj case ne treba vec moze sve do 4-tog reda ( ili kasnije osmog # ) 01424 switch(this->rank()) 01425 { 01426 case 0: 01427 { 01428 // ::printf(" trace=%.4e ", val(1)); 01429 // ::printf("\n"); 01430 tr = cval(1); 01431 break; 01432 } 01433 01434 case 1: 01435 { 01436 if(dim()[0] != 1) 01437 { 01438 ::printf("\a\nERROR in trace function : not a squared 1-st rank BJtensor\n"); 01439 ::exit( 1 ); 01440 } 01441 tr = cval(1); 01442 break; 01443 } 01444 01445 case 2: 01446 { 01447 if(dim()[0] != dim()[1]) 01448 { 01449 ::printf("\a\nERROR in trace function : not a sqared 2nd-rank BJtensor\n"); 01450 ::exit( 1 ); 01451 } 01452 for ( int i2=1 ; i2<=dim()[0] ; i2++ ) 01453 tr += cval(i2, i2); 01454 break; 01455 } 01456 01457 case 3: 01458 { 01459 if( dim()[0] != dim()[1] || 01460 dim()[1] != dim()[2] || 01461 dim()[2] != dim()[0] ) 01462 { 01463 ::printf("\a\nERROR in trace function : not a sqared 3nd-rank BJtensor\n"); 01464 ::exit( 1 ); 01465 } 01466 for ( int i3=1 ; i3<=dim()[0] ; i3++ ) 01467 tr += cval(i3, i3, i3); 01468 break; 01469 } 01470 01471 case 4: 01472 { 01473 if( dim()[0] != dim()[1] || 01474 dim()[1] != dim()[2] || 01475 dim()[2] != dim()[3] || 01476 dim()[3] != dim()[0] ) 01477 { 01478 ::printf("\a\nERROR in trace function : not a sqared 4nd-rank BJtensor\n"); 01479 ::exit( 1 ); 01480 } 01481 for ( int i4=1 ; i4<=dim()[0] ; i4++ ) 01482 tr += cval(i4, i4, i4, i4); 01483 break; 01484 } 01485 } 01486 01487 //.......... null_indices(); 01488 01489 return tr; 01490 } 01491 01492 01493 01494 //-----//############################################################################## 01495 //-----// scalar multiplication (nDarray * scalar) 01496 //-----nDarray nDarray::operator*( double rval) 01497 //----- { 01498 //-----// construct nDarray using the same control numbers as for the 01499 //-----// original one. 01500 //----- nDarray mult(this->pc_nDarray_rep->nDarray_rank, this->pc_nDarray_rep->dim, 0.0); 01501 //----- 01502 //----- switch(pc_nDarray_rep->nDarray_rank) 01503 //----- { 01504 //----- case 0: 01505 //----- { 01506 //----- mult.val(1) = val(1) * rval; 01507 //----- break; 01508 //----- } 01509 //----- 01510 //----- case 1: 01511 //----- { 01512 //----- for ( int i1=1 ; i1<=this->pc_nDarray_rep->dim[0] ; i1++ ) 01513 //----- mult.val(i1) = val(i1) * rval; 01514 //----- break; 01515 //----- } 01516 //----- 01517 //----- case 2: 01518 //----- { 01519 //----- for ( int i2=1 ; i2<=this->pc_nDarray_rep->dim[0] ; i2++ ) 01520 //----- for ( int j2=1 ; j2<=this->pc_nDarray_rep->dim[1] ; j2++ ) 01521 //----- mult.val(i2, j2) = val(i2, j2) * rval; 01522 //----- break; 01523 //----- } 01524 //----- 01525 //----- case 3: 01526 //----- { 01527 //----- for ( int i3=1 ; i3<=this->pc_nDarray_rep->dim[0] ; i3++ ) 01528 //----- for ( int j3=1 ; j3<=this->pc_nDarray_rep->dim[1] ; j3++ ) 01529 //----- for ( int k3=1 ; k3<=this->pc_nDarray_rep->dim[2] ; k3++ ) 01530 //----- mult.val(i3, j3, k3) = val(i3, j3, k3) * rval; 01531 //----- break; 01532 //----- } 01533 //----- 01534 //----- case 4: 01535 //----- { 01536 //----- for ( int i4=1 ; i4<=this->pc_nDarray_rep->dim[0] ; i4++ ) 01537 //----- for ( int j4=1 ; j4<=this->pc_nDarray_rep->dim[1] ; j4++ ) 01538 //----- for ( int k4=1 ; k4<=this->pc_nDarray_rep->dim[2] ; k4++ ) 01539 //----- for ( int l4=1 ; l4<=this->pc_nDarray_rep->dim[3] ; l4++ ) 01540 //----- mult.val(i4,j4,k4,l4)=val(i4,j4,k4,l4)*rval; 01541 //----- break; 01542 //----- } 01543 //----- } 01544 //----- 01545 //----- return mult; 01546 //----- } 01547 //----- 01548 01549 //############################################################################## 01550 // nDarray comparisson // nDarray comparisson 01551 int nDarray::operator==( nDarray & rval) // returns 1 if they are same 01552 { // returns 0 if they are not 01553 int true_or_not = 1; // suppose that they are the same 01554 01555 double sqrt_d_macheps = sqrt(d_macheps()); 01556 // double qub_d_macheps = pow(d_macheps(),(1.0/3.0)); 01557 01558 double tolerance = sqrt_d_macheps; 01559 01560 int this_rank_of_nDarray = this->pc_nDarray_rep->nDarray_rank; 01561 int rval_rank_of_nDarray = rval.pc_nDarray_rep->nDarray_rank; 01562 01563 if(this_rank_of_nDarray != rval_rank_of_nDarray) 01564 { 01565 ::printf("\a\nnDarrays of different ranks: comparisson not possible\n"); 01566 ::exit ( 1 ); 01567 } 01568 01569 for ( int i=0 ; i<this_rank_of_nDarray ; i++ ) 01570 if (this->pc_nDarray_rep->dim[i] != rval.pc_nDarray_rep->dim[i] ) 01571 { 01572 ::fprintf(stderr,"\a\nDimension discrepancy in operator==\n\ 01573 this->pc_nDarray_rep->dim[%d]=%d\n\ 01574 arg.pc_nDarray_rep->dim[%d]=%d\n", 01575 i,this->pc_nDarray_rep->dim[i], 01576 i,rval.pc_nDarray_rep->dim[i]); 01577 ::exit(1); 01578 } 01579 //..// construct nDarray using the same control numbers as for the 01580 //..// original one . 01581 //.. nDarray add(pc_nDarray_rep->nDarray_rank, pc_nDarray_rep->dim, 0.0); 01582 01583 switch(pc_nDarray_rep->nDarray_rank) 01584 { 01585 case 0: 01586 { 01587 if ( fabs( val(1)-rval.val(1) ) >= tolerance) 01588 true_or_not = 0; 01589 break; 01590 } 01591 01592 case 1: 01593 { 01594 for ( int i1=1 ; i1<=this->pc_nDarray_rep->dim[0] ; i1++ ) 01595 if ( fabs( val(i1)-rval.val(i1) ) >= tolerance) 01596 true_or_not = 0; 01597 break; 01598 } 01599 01600 case 2: 01601 { 01602 for ( int i2=1 ; i2<=this->pc_nDarray_rep->dim[0] ; i2++ ) 01603 for ( int j2=1 ; j2<=this->pc_nDarray_rep->dim[1] ; j2++ ) 01604 if ( fabs( val(i2,j2)-rval.val(i2,j2) ) >= tolerance) 01605 true_or_not = 0; 01606 break; 01607 } 01608 01609 case 3: 01610 { 01611 for ( int i3=1 ; i3<=this->pc_nDarray_rep->dim[0] ; i3++ ) 01612 for ( int j3=1 ; j3<=this->pc_nDarray_rep->dim[1] ; j3++ ) 01613 for ( int k3=1 ; k3<=this->pc_nDarray_rep->dim[2] ; k3++ ) 01614 if (fabs(val(i3,j3,k3)-rval.val(i3,j3,k3))>=tolerance) 01615 true_or_not = 0; 01616 break; 01617 } 01618 01619 case 4: 01620 { 01621 for ( int i4=1 ; i4<=this->pc_nDarray_rep->dim[0] ; i4++ ) 01622 for ( int j4=1 ; j4<=this->pc_nDarray_rep->dim[1] ; j4++ ) 01623 for ( int k4=1 ; k4<=this->pc_nDarray_rep->dim[2] ; k4++ ) 01624 for ( int l4=1 ; l4<=this->pc_nDarray_rep->dim[3] ; l4++ ) 01625 if(fabs(val(i4,j4,k4,l4)-rval.val(i4,j4,k4,l4))>=tolerance) 01626 true_or_not = 0; 01627 break; 01628 } 01629 } 01630 01631 return true_or_not; 01632 } 01633 01634 01635 01636 //############################################################################## 01637 // nDarray print function 01638 void nDarray::print(char *name , char *msg) const 01639 { 01640 if (*msg) ::printf("%s\n",msg); 01641 01642 switch(pc_nDarray_rep->nDarray_rank) 01643 { 01644 case 0: 01645 { 01646 ::printf("%s(1)=%+8.4e ", name, cval(1)); 01647 ::printf("\n"); 01648 break; 01649 } 01650 01651 case 1: 01652 { 01653 for ( int i=1 ; i<=pc_nDarray_rep->dim[0]; i++ ) 01654 { 01655 ::printf("%s(%2d)=%+8.4e ", name, i, cval(i)); 01656 ::printf("\n"); 01657 } 01658 break; 01659 } 01660 01661 case 2: 01662 { 01663 for ( int i2=1 ; i2<=pc_nDarray_rep->dim[0] ; i2++ ) 01664 { 01665 for ( int j2=1 ; j2<=pc_nDarray_rep->dim[1] ; j2++ ) 01666 { 01667 ::printf("%s(%2d,%2d)=%+12.8e ", name, i2, j2, cval(i2, j2)); 01668 } 01669 ::printf("\n"); 01670 } 01671 break; 01672 } 01673 01674 case 3: 01675 { 01676 for ( int i3=1 ; i3<=pc_nDarray_rep->dim[0] ; i3++ ) 01677 for ( int j3=1 ; j3<=pc_nDarray_rep->dim[1] ; j3++ ) 01678 { 01679 for ( int k3=1 ; k3<=pc_nDarray_rep->dim[2] ; k3++ ) 01680 { 01681 ::printf("%s(%d,%d,%d)=%+8.4e ", name, i3, j3, k3, 01682 cval(i3, j3, k3)); 01683 } 01684 ::printf("\n"); 01685 } 01686 break; 01687 } 01688 01689 case 4: 01690 { 01691 for ( int i4=1 ; i4<=pc_nDarray_rep->dim[0] ; i4++ ) 01692 for ( int j4=1 ; j4<=pc_nDarray_rep->dim[1] ; j4++ ) 01693 for ( int k4=1 ; k4<=pc_nDarray_rep->dim[2] ; k4++ ) 01694 { 01695 for ( int l4=1 ; l4<=pc_nDarray_rep->dim[3] ; l4++ ) 01696 { 01697 ::printf("%s(%d,%d,%d,%d)=%+8.4e ",name,i4,j4,k4,l4, 01698 cval(i4, j4, k4, l4)); 01699 } 01700 ::printf("\n"); 01701 } 01702 break; 01703 } 01704 } 01706 // 01707 } 01708 01709 //############################################################################## 01710 // nDarray print function 01711 void nDarray::printshort(char *msg) const 01712 { 01713 if (*msg) ::printf("%s\n",msg); 01714 01715 switch(pc_nDarray_rep->nDarray_rank) 01716 { 01717 case 0: 01718 { 01719 ::printf("%+6.2e ",cval(1)); 01720 ::printf("\n"); 01721 break; 01722 } 01723 01724 case 1: 01725 { 01726 for ( int i=1 ; i<=pc_nDarray_rep->dim[0]; i++ ) 01727 { 01728 ::printf("%+6.2e ",cval(i)); 01729 ::printf("\n"); 01730 } 01731 break; 01732 } 01733 01734 case 2: 01735 { 01736 for ( int i2=1 ; i2<=pc_nDarray_rep->dim[0] ; i2++ ) 01737 { 01738 for ( int j2=1 ; j2<=pc_nDarray_rep->dim[1] ; j2++ ) 01739 { 01740 ::printf("%+6.2e ", cval(i2, j2)); 01741 } 01742 ::printf("\n"); 01743 } 01744 break; 01745 } 01746 01747 case 3: 01748 { 01749 for ( int i3=1 ; i3<=pc_nDarray_rep->dim[0] ; i3++ ) 01750 for ( int j3=1 ; j3<=pc_nDarray_rep->dim[1] ; j3++ ) 01751 { 01752 for ( int k3=1 ; k3<=pc_nDarray_rep->dim[2] ; k3++ ) 01753 { 01754 ::printf("%+6.2e ", cval(i3, j3, k3)); 01755 } 01756 ::printf("\n"); 01757 } 01758 break; 01759 } 01760 01761 case 4: 01762 { 01763 for ( int i4=1 ; i4<=pc_nDarray_rep->dim[0] ; i4++ ) 01764 for ( int j4=1 ; j4<=pc_nDarray_rep->dim[1] ; j4++ ) 01765 for ( int k4=1 ; k4<=pc_nDarray_rep->dim[2] ; k4++ ) 01766 { 01767 for ( int l4=1 ; l4<=pc_nDarray_rep->dim[3] ; l4++ ) 01768 { 01769 ::printf("%+6.2e ", cval(i4, j4, k4, l4)); 01770 } 01771 ::printf("\n"); 01772 } 01773 break; 01774 } 01775 } 01777 // 01778 } 01779 01780 //############################################################################## 01781 // nDarray print function for mathematica 01782 void nDarray::mathprint(void) const 01783 { 01784 01785 switch(pc_nDarray_rep->nDarray_rank) 01786 { 01787 case 0: 01788 { 01789 ::printf("{"); 01790 ::printf("%12f, ", cval(1)); 01791 ::printf("},\n"); 01792 break; 01793 } 01794 01795 case 1: 01796 { 01797 ::printf("{"); 01798 int i = 1; 01799 for ( i=1 ; i<=pc_nDarray_rep->dim[0]; i++ ) 01800 { 01801 ::printf("%12f, ", cval(i)); 01802 ::printf("\n"); 01803 if (i<pc_nDarray_rep->dim[0]) ::printf(", "); 01804 if (i==pc_nDarray_rep->dim[0]) ::printf(" \n"); 01805 } 01806 if (i<(pc_nDarray_rep->dim[0]+1)) ::printf("},\n"); 01807 if (i==(pc_nDarray_rep->dim[0]+1)) ::printf("}\n"); 01808 break; 01809 } 01810 01811 case 2: 01812 { 01813 ::printf("{\n"); 01814 for ( int i2=1 ; i2<=pc_nDarray_rep->dim[0] ; i2++ ) 01815 { 01816 if (pc_nDarray_rep->dim[1]!=1) 01817 { 01818 ::printf("{"); 01819 } 01820 for ( int j2=1 ; j2<=pc_nDarray_rep->dim[1] ; j2++ ) 01821 { 01822 ::printf("%12f", cval(i2, j2)); 01823 if (j2<pc_nDarray_rep->dim[1] ) 01824 { 01825 ::printf(", "); 01826 } 01827 // if (j2==pc_nDarray_rep->dim[1]) ::printf(""); 01828 } 01829 // if (i2<pc_nDarray_rep->dim[1] && pc_nDarray_rep->dim[1]!=1 ) 01830 // { 01831 // ::printf(",\n"); 01832 // } 01833 if ( pc_nDarray_rep->dim[1]==1 && i2<pc_nDarray_rep->dim[0] ) 01834 { 01835 ::printf(", "); 01836 } 01837 if (i2<pc_nDarray_rep->dim[1] ) 01838 { 01839 ::printf("},\n"); 01840 } 01841 if (i2==pc_nDarray_rep->dim[1] && pc_nDarray_rep->dim[1]!=1 ) 01842 { 01843 ::printf("}\n"); 01844 } 01845 } 01846 // if (i2<(pc_nDarray_rep->dim[0]+1)) ::printf("},\n"); 01847 // if (i2==(pc_nDarray_rep->dim[0]+1)) 01848 // { 01849 ::printf("}\n"); 01850 // } 01851 break; 01852 } 01853 01854 case 3: 01855 { 01856 for ( int i3=1 ; i3<=pc_nDarray_rep->dim[0] ; i3++ ) 01857 for ( int j3=1 ; j3<=pc_nDarray_rep->dim[1] ; j3++ ) 01858 { 01859 for ( int k3=1 ; k3<=pc_nDarray_rep->dim[2] ; k3++ ) 01860 { 01861 ::printf("%12f, ", cval(i3, j3, k3)); 01862 } 01863 ::printf("\n"); 01864 } 01865 break; 01866 } 01867 01868 case 4: 01869 { 01870 for ( int i4=1 ; i4<=pc_nDarray_rep->dim[0] ; i4++ ) 01871 for ( int j4=1 ; j4<=pc_nDarray_rep->dim[1] ; j4++ ) 01872 for ( int k4=1 ; k4<=pc_nDarray_rep->dim[2] ; k4++ ) 01873 { 01874 for ( int l4=1 ; l4<=pc_nDarray_rep->dim[3] ; l4++ ) 01875 { 01876 ::printf("%12f, ", cval(i4, j4, k4, l4)); 01877 } 01878 ::printf("\n"); 01879 } 01880 break; 01881 } 01882 } 01884 // 01885 } 01886 01887 01888 //############################################################################## 01889 // nDarray Frobenius_norm function // return the Frobenius norm of 01890 double nDarray::Frobenius_norm(void) // BJmatrix, BJtensor, BJvector 01891 { // see Dennis & Schnabel page 43 01892 double FN = 0.0; 01893 01894 switch(pc_nDarray_rep->nDarray_rank) 01895 { 01896 case 0: 01897 { 01898 FN = (val(1))*(val(1)); 01899 break; 01900 } 01901 01902 case 1: 01903 { 01904 for ( int i=1 ; i<=pc_nDarray_rep->dim[0]; i++ ) 01905 { 01906 FN = FN + (val(i))*(val(i)); 01907 } 01908 break; 01909 } 01910 01911 case 2: 01912 { 01913 for ( int i2=1 ; i2<=pc_nDarray_rep->dim[0] ; i2++ ) 01914 { 01915 for ( int j2=1 ; j2<=pc_nDarray_rep->dim[1] ; j2++ ) 01916 { 01917 FN = FN + (val(i2,j2))*(val(i2,j2)); 01918 } 01919 } 01920 break; 01921 } 01922 01923 case 3: 01924 { 01925 for ( int i3=1 ; i3<=pc_nDarray_rep->dim[0] ; i3++ ) 01926 for ( int j3=1 ; j3<=pc_nDarray_rep->dim[1] ; j3++ ) 01927 { 01928 for ( int k3=1 ; k3<=pc_nDarray_rep->dim[2] ; k3++ ) 01929 { 01930 FN = FN + (val(i3,j3,k3))*(val(i3,j3,k3)); 01931 } 01932 } 01933 break; 01934 } 01935 01936 case 4: 01937 { 01938 for ( int i4=1 ; i4<=pc_nDarray_rep->dim[0] ; i4++ ) 01939 for ( int j4=1 ; j4<=pc_nDarray_rep->dim[1] ; j4++ ) 01940 for ( int k4=1 ; k4<=pc_nDarray_rep->dim[2] ; k4++ ) 01941 { 01942 for ( int l4=1 ; l4<=pc_nDarray_rep->dim[3] ; l4++ ) 01943 { 01944 FN = FN + (val(i4, j4, k4, l4))*(val(i4, j4, k4, l4)); 01945 } 01946 } 01947 break; 01948 } 01949 } 01950 01951 return sqrt(FN); 01952 01953 } 01954 01955 //############################################################################## 01956 // nDarray General_norm function // return the General p-th norm of 01957 double nDarray::General_norm(double p) // BJmatrix, BJtensor, BJvector 01958 { // see Dennis & Schnabel page 42 01959 double N = 0.0; 01960 01961 switch(pc_nDarray_rep->nDarray_rank) 01962 { 01963 case 0: 01964 { 01965 N = pow(fabs(val(1)),p); 01966 break; 01967 } 01968 01969 case 1: 01970 { 01971 for ( int i=1 ; i<=pc_nDarray_rep->dim[0]; i++ ) 01972 { 01973 N = N + pow(fabs(val(i)),p); 01974 } 01975 break; 01976 } 01977 01978 case 2: 01979 { 01980 for ( int i2=1 ; i2<=pc_nDarray_rep->dim[0] ; i2++ ) 01981 { 01982 for ( int j2=1 ; j2<=pc_nDarray_rep->dim[1] ; j2++ ) 01983 { 01984 N = N + pow(fabs(val(i2,j2)),p); 01985 } 01986 } 01987 break; 01988 } 01989 01990 case 3: 01991 { 01992 for ( int i3=1 ; i3<=pc_nDarray_rep->dim[0] ; i3++ ) 01993 for ( int j3=1 ; j3<=pc_nDarray_rep->dim[1] ; j3++ ) 01994 { 01995 for ( int k3=1 ; k3<=pc_nDarray_rep->dim[2] ; k3++ ) 01996 { 01997 N = N + pow(fabs(val(i3,j3,k3)),p); 01998 } 01999 } 02000 break; 02001 } 02002 02003 case 4: 02004 { 02005 for ( int i4=1 ; i4<=pc_nDarray_rep->dim[0] ; i4++ ) 02006 for ( int j4=1 ; j4<=pc_nDarray_rep->dim[1] ; j4++ ) 02007 for ( int k4=1 ; k4<=pc_nDarray_rep->dim[2] ; k4++ ) 02008 { 02009 for ( int l4=1 ; l4<=pc_nDarray_rep->dim[3] ; l4++ ) 02010 { 02011 N = N + pow(fabs(val(i4,j4,k4,l4)),p); 02012 } 02013 } 02014 break; 02015 } 02016 } 02017 02018 return pow(N,(1.0/p)); 02019 02020 } 02021 02022 02023 02024 //############################################################################## 02025 int nDarray::number_of_zeros() const // number of members that are 02026 { // smaller than sqrt(macheps) 02027 int n = 0; 02028 double machepsilon = d_macheps(); 02029 double tolerance = sqrt(machepsilon); 02030 02031 for (int j=0 ; j<this->pc_nDarray_rep->total_numb ; j++) 02032 if ( this->pc_nDarray_rep->pd_nDdata[j] <= tolerance ) 02033 { 02034 n++; 02035 } 02036 02037 return n; 02038 } 02039 02040 02041 // prebacen u nDarray 14 oktobra 1996 02042 //############################################################################# 02043 nDarray nDarray::eigenvalues(void) 02044 { 02045 int rows = this->dim(1); 02046 int cols = this->dim(2); 02047 if ( rows != cols && this->rank() != 2 ) 02048 { 02049 ::printf("rows!=cols in eigenvalues and rank != 2 \n"); 02050 ::exit(1); 02051 } 02052 02053 // BJvector EV(rows, 0.0); 02054 const int pdim[] = {rows}; 02055 nDarray EV(1, pdim, 0.0); 02056 02057 // najbezbolinije da stvarno napravim dvodimenzioni niz pa da ga kopiram u 'a' 02058 // PAZI oni u NRC rade kao u FORTRANU dakle nizovi od 1 do n 02059 double ** a = new double *[rows+1]; 02060 if ( !a ) {::printf("memory exhausted for **a \n"); ::exit(1);} 02061 for ( int i=0 ; i<(rows+1) ; i++ ) 02062 { 02063 a[i] = new double[rows+1]; 02064 if ( !a[i] ) {::printf("memory exhausted for *a \n"); ::exit(1);} 02065 } 02066 double * d = new double [rows+1]; 02067 double * e = new double [rows+1]; 02068 02069 for ( int j=0 ; j<rows ; j++) 02070 for ( int k=0 ; k<rows ; k++) 02071 { 02072 a[j+1][k+1] = this->cval(j+1,k+1); 02073 } 02074 02075 02076 // Householder reduction of a real symmetric BJmatrix see NRC page 373 02077 tred2( a, rows, d, e); 02078 // QL algorithm with implicit shift see NRC page 380 02079 tqli( d, e, rows, a); 02080 // sort eigenvalues see NRC page 366 02081 eigsrt( d, a, rows); 02082 02083 for ( int l=0 ; l<rows ; l++ ) 02084 { 02085 EV.val(l+1) = d[l+1]; 02086 } 02087 02088 // delete new - ed arrays 02089 // for ( int i1=0 ; i1<rows ; i1++ ) 02090 // { 02091 // delete a[i1]; 02092 // } 02093 for ( int ii=0 ; ii<(rows+1) ; ii++ ) 02094 { 02095 delete [] a[ii]; 02096 } 02097 delete [] a; 02098 delete [] d; 02099 delete [] e; 02100 02101 return EV; 02102 } 02103 02104 //############################################################################# 02105 nDarray nDarray::eigenvectors(void) 02106 { 02107 int rows = this->dim(1); 02108 int cols = this->dim(2); 02109 if ( rows != cols && this->rank() != 2 ) 02110 { 02111 ::printf("rows!=cols in eigenBJvectors and rank != 2 \n"); 02112 ::exit(1); 02113 } 02114 02115 // BJmatrix EV(rows, rows, 0.0); 02116 const int pdim[] = {rows, rows}; 02117 nDarray EV(2, pdim, 0.0); 02118 // BJmatrix temp( rows, rows, rows, this->data() ); 02119 02120 // najbezbolinije da stvarno napravim dvodimenzioni niz pa da ga kopiram u 'a' 02121 // PAZI oni u NRC rade kao u FORTRANU dakle nizovi od 1 - n 02122 double ** a = new double *[rows+1]; 02123 if ( !a ) {::printf("memory exhausted for **a \n"); ::exit(1);} 02124 for ( int i=0 ; i<(rows+1) ; i++ ) 02125 { 02126 a[i] = new double[rows+1]; 02127 if ( !a[i] ) {::printf("memory exhausted for *a \n"); ::exit(1);} 02128 } 02129 double * d = new double [rows+1]; 02130 double * e = new double [rows+1]; 02131 02132 for ( int j=0 ; j<rows ; j++) 02133 for ( int k=0 ; k<rows ; k++) 02134 { 02135 a[j+1][k+1] = this->cval(j+1,k+1); 02136 } 02137 // Householder reduction of a real symmetric BJmatrix see NRC page 373 02138 tred2( a, rows, d, e); 02139 // QL algorithm with implicit shift see NRC page 380 02140 tqli( d, e, rows, a); 02141 // sort eigenvalues see NRC page 366 02142 eigsrt( d, a, rows); 02143 02144 for ( int l=0 ; l<rows ; l++ ) 02145 for ( int l1=0; l1<rows ; l1++ ) 02146 { 02147 EV.val(l+1,l1+1) = a[l+1][l1+1]; 02148 } 02149 02150 // delete new - ed arrays 02151 for ( int i1=0 ; i1<(rows+1) ; i1++ ) 02152 { 02153 delete [] a[i1]; 02154 } 02155 delete [] a; 02156 delete [] d; 02157 delete [] e; 02158 02159 return EV; 02160 } 02161 02162 02163 //############################################################################## 02164 nDarray nDarray::nDsqrt(void) // returns all elements square root 02165 { 02166 nDarray newnD( this->rank(), this->pc_nDarray_rep->dim, this->data()); 02167 02168 for ( int i=0 ; i<pc_nDarray_rep->total_numb ; i++ ) 02169 newnD.pc_nDarray_rep->pd_nDdata[i] = sqrt(this->pc_nDarray_rep->pd_nDdata[i]) ; 02170 02171 return newnD; 02172 02173 } 02174 02175 02176 02177 02178 // prebacen u nDarray 14 oktobra 1996 02179 //############################################################################# 02180 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a)) 02181 02182 void nDarray::tqli(double * d, double * e, int n, double ** z) 02183 { 02184 int m,l,iter,i,k; 02185 double s,r,p,g,f,dd,c,b; 02186 // void nrerror(); 02187 02188 for (i=2;i<=n;i++) e[i-1]=e[i]; 02189 e[n]=0.0; 02190 for (l=1;l<=n;l++) { 02191 iter=0; 02192 do { 02193 for (m=l;m<=n-1;m++) { 02194 dd=fabs(d[m])+fabs(d[m+1]); 02195 if (fabs(e[m])+dd == dd) break; 02196 } 02197 if (m != l) { 02198 if (iter++ == 30) { ::printf("Too many iterations in TQLI\n"); ::exit(1); } 02199 g=(d[l+1]-d[l])/(2.0*e[l]); 02200 r=sqrt((g*g)+1.0); 02201 g=d[m]-d[l]+e[l]/(g+SIGN(r,g)); 02202 s=c=1.0; 02203 p=0.0; 02204 for (i=m-1;i>=l;i--) { 02205 f=s*e[i]; 02206 b=c*e[i]; 02207 if (fabs(f) >= fabs(g)) { 02208 c=g/f; 02209 r=sqrt((c*c)+1.0); 02210 e[i+1]=f*r; 02211 c *= (s=1.0/r); 02212 } else { 02213 s=f/g; 02214 r=sqrt((s*s)+1.0); 02215 e[i+1]=g*r; 02216 s *= (c=1.0/r); 02217 } 02218 g=d[i+1]-p; 02219 r=(d[i]-g)*s+2.0*c*b; 02220 p=s*r; 02221 d[i+1]=g+p; 02222 g=c*r-b; 02223 /* Next loop can be omitted if eigenBJvectors not wanted */ 02224 for (k=1;k<=n;k++) { 02225 f=z[k][i+1]; 02226 z[k][i+1]=s*z[k][i]+c*f; 02227 z[k][i]=c*z[k][i]-s*f; 02228 } 02229 } 02230 d[l]=d[l]-p; 02231 e[l]=g; 02232 e[m]=0.0; 02233 } 02234 } while (m != l); 02235 } 02236 } 02237 02238 02239 //############################################################################# 02240 void nDarray::tred2(double ** a, int n, double * d, double * e) 02241 { 02242 int l,k,j,i; 02243 double scale,hh,h,g,f; 02244 02245 for (i=n;i>=2;i--) { 02246 l=i-1; 02247 h=scale=0.0; 02248 if (l > 1) { 02249 for (k=1;k<=l;k++) 02250 scale += fabs(a[i][k]); 02251 if (scale == 0.0) 02252 e[i]=a[i][l]; 02253 else { 02254 for (k=1;k<=l;k++) { 02255 a[i][k] /= scale; 02256 h += a[i][k]*a[i][k]; 02257 } 02258 f=a[i][l]; 02259 g = f>0 ? -sqrt(h) : sqrt(h); 02260 e[i]=scale*g; 02261 h -= f*g; 02262 a[i][l]=f-g; 02263 f=0.0; 02264 for (j=1;j<=l;j++) { 02265 /* Next statement can be omitted if eigenBJvectors not wanted */ 02266 a[j][i]=a[i][j]/h; 02267 g=0.0; 02268 for (k=1;k<=j;k++) 02269 g += a[j][k]*a[i][k]; 02270 for (k=j+1;k<=l;k++) 02271 g += a[k][j]*a[i][k]; 02272 e[j]=g/h; 02273 f += e[j]*a[i][j]; 02274 } 02275 hh=f/(h+h); 02276 for (j=1;j<=l;j++) { 02277 f=a[i][j]; 02278 e[j]=g=e[j]-hh*f; 02279 for (k=1;k<=j;k++) 02280 a[j][k] -= (f*e[k]+g*a[i][k]); 02281 } 02282 } 02283 } else 02284 e[i]=a[i][l]; 02285 d[i]=h; 02286 } 02287 /* Next statement can be omitted if eigenBJvectors not wanted */ 02288 d[1]=0.0; 02289 e[1]=0.0; 02290 /* Contents of this loop can be omitted if eigenBJvectors not 02291 wanted except for statement d[i]=a[i][i]; */ 02292 for (i=1;i<=n;i++) { 02293 l=i-1; 02294 if (d[i]) { 02295 for (j=1;j<=l;j++) { 02296 g=0.0; 02297 for (k=1;k<=l;k++) 02298 g += a[i][k]*a[k][j]; 02299 for (k=1;k<=l;k++) 02300 a[k][j] -= g*a[k][i]; 02301 } 02302 } 02303 d[i]=a[i][i]; 02304 a[i][i]=1.0; 02305 for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0; 02306 } 02307 } 02308 02309 02310 02311 //############################################################################# 02312 void nDarray::eigsrt(double * d, double ** v, int n) 02313 { 02314 int k,j,i; 02315 double p; 02316 02317 for (i=1;i<n;i++) { 02318 p=d[k=i]; 02319 for (j=i+1;j<=n;j++) 02320 if (d[j] >= p) p=d[k=j]; 02321 if (k != i) { 02322 d[k]=d[i]; 02323 d[i]=p; 02324 for (j=1;j<=n;j++) { 02325 p=v[j][i]; 02326 v[j][i]=v[j][k]; 02327 v[j][k]=p; 02328 } 02329 } 02330 } 02331 } 02332 02333 02334 02335 02336 02337 //############################################################################## 02338 // very private part 02339 //############################################################################## 02340 double * nDarray::data(void) const 02341 { 02342 return this->pc_nDarray_rep->pd_nDdata; 02343 } 02344 02345 void nDarray::set_data_pointer(double * data_pointer) 02346 { 02347 this->pc_nDarray_rep->pd_nDdata = data_pointer; 02348 } 02349 02350 int nDarray::rank(void) const 02351 { 02352 return this->pc_nDarray_rep->nDarray_rank; 02353 } 02354 02355 void nDarray::rank(int nDarray_rank) 02356 { 02357 this->pc_nDarray_rep->nDarray_rank = nDarray_rank; 02358 } 02359 02360 long int nDarray::total_number(void) const 02361 { 02362 return this->pc_nDarray_rep->total_numb; 02363 } 02364 02365 void nDarray::total_number(long int number) 02366 { 02367 this->pc_nDarray_rep->total_numb = number; 02368 } 02369 02370 int * nDarray::dim(void) const 02371 { 02372 return this->pc_nDarray_rep->dim; 02373 } 02374 02375 int & nDarray::get_dim_pointer(void) const 02376 { 02377 return this->pc_nDarray_rep->dim[0]; 02378 } 02379 02380 void nDarray::set_dim_pointer(int * dim_pointer) 02381 { 02382 this->pc_nDarray_rep->dim = dim_pointer; 02383 } 02384 02385 int nDarray::dim(int which) const 02386 { 02387 return this->pc_nDarray_rep->dim[which-1]; 02388 } 02389 02390 int nDarray::reference_count(int up_down) 02391 { 02392 this->pc_nDarray_rep->n += up_down; 02393 return(this->pc_nDarray_rep->n); 02394 } 02395 02396 void nDarray::set_reference_count(int ref_count) 02397 { 02398 this->pc_nDarray_rep->n=ref_count; 02399 } 02400 02401 02402 //tempOUT//############################################################################## 02403 // TENSOR_REP_CC 02404 // ####################### 02405 // memory manager part 02406 // overloading operator new in nDarray::nDarray_rep class ################## 02407 void * nDarray_rep::operator new(size_t s) 02408 { // see C++ reference manual by 02409 void *void_pointer; // ELLIS and STROUSTRUP page 283. 02410 void_pointer = ::operator new(s); // and ECKEL page 529. 02411 // ::printf("\nnew pointer %p of size %d\n",void_pointer,s); 02412 if (!void_pointer) 02413 { 02414 ::fprintf(stderr,"\a\nInsufficient memory\n"); 02415 ::exit(1); 02416 } 02417 return void_pointer; 02418 } 02419 02420 // overloading operator delete in nDarray::nDarray_rep class ################## 02421 void nDarray_rep::operator delete(void *p) 02422 { // see C++ reference manual by 02423 // ELLIS and STROUSTRUP page 283. 02424 // and ECKEL page 529. 02425 // ::printf("deleted pointer %p\n",p); 02426 ::operator delete(p); 02427 } 02428 02429 02430 02431 02432 #endif 02433 02434
Generated on Mon Oct 23 15:05:24 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents