00001
00002
00003
00004
00005
00006
00007 #include "AdmittanceNetwork_PAR1.h"
00008 #include <iostream>
00009 #include <cassert>
00010
00011 #include <mpi.h>
00012 #include <zmumps_c.h>
00013 #ifndef ZMUMPS_COMPLEX
00014 #define ZMUMPS_COMPLEX mumps_double_complex
00015 #endif
00016 #include <map>
00017 using namespace std;
00018
00019
00020 #define TESTMAIN testmain
00021 #include "GreedyBW.cpp"
00022
00023
00024 typedef map<int,Complex> SparseComplexVector;
00025 typedef map<int, SparseComplexVector> SparseComplexMatrix;
00026
00027
00028 struct PAR1Data
00029 {
00030 PAR1Data() : id(), inited(false), ncmds(0),
00031 nsolves(0), nfactors(0),
00032 analysistime(0), solvetime(0), factortime(0) {}
00033 ZMUMPS_STRUC_C id;
00034 bool inited;
00035 int ncmds, nsolves, nfactors;
00036 double analysistime, solvetime, factortime;
00037 };
00038
00039
00040 class AdmittanceNetwork_PAR1::Data
00041 {
00042 public:
00043 Data(int n) : N(n), Y(), perm() {}
00044 virtual ~Data() {}
00045
00046 public:
00047 const int N;
00048 SparseComplexMatrix Y;
00049 Labels perm;
00050
00051 struct MPIData
00052 {
00053 int rank, nranks;
00054 MPI_Comm comm;
00055 } mpi;
00056
00057 PAR1Data par1;
00058 };
00059
00060
00061 AdmittanceNetwork_PAR1::AdmittanceNetwork_PAR1(int N):
00062 dirty(true)
00063 {
00064 data = new Data( N );
00065 data->mpi.comm = MPI_COMM_WORLD;
00066 MPI_Comm_rank( data->mpi.comm, &data->mpi.rank );
00067 MPI_Comm_size( data->mpi.comm, &data->mpi.nranks );
00068 initlabels( data->perm, N, true );
00069 }
00070
00071
00072 AdmittanceNetwork_PAR1::~AdmittanceNetwork_PAR1()
00073 {
00074 delete data;
00075 }
00076
00077
00078 void AdmittanceNetwork_PAR1::set(int i, int j, Complex y)
00079 {
00080 dirty = true;
00081 data->Y[i][j] = y;
00082 }
00083
00084
00085 Complex AdmittanceNetwork_PAR1::get(int i, int j) const
00086 {
00087 Complex val(0,0);
00088 SparseComplexMatrix *mat = &data->Y;
00089 SparseComplexMatrix::const_iterator mit = mat->find( i );
00090 if( mit != mat->end() )
00091 {
00092 const SparseComplexVector &vec = mit->second;
00093 SparseComplexVector::const_iterator vit = vec.find( j );
00094 if( vit != vec.end() )
00095 {
00096 val = vit->second;
00097 }
00098 }
00099 return val;
00100 }
00101
00102
00103 void AdmittanceNetwork_PAR1::add_line(int i, int j, Complex y)
00104 {
00105 set(i,i,get(i,i)+y);
00106 set(j,j,get(j,j)+y);
00107 set(i,j,get(i,j)-y);
00108 set(j,i,get(j,i)-y);
00109 }
00110
00111
00112 void AdmittanceNetwork_PAR1::add_line(int i, int j, Complex y, double turns,
00113 double phase_shift, int side)
00114 {
00115 Complex a = polar(turns,phase_shift);
00116 set(side,side,get(side,side)+a*a*y);
00117 int off_side = i;
00118 if (side == off_side) off_side = j;
00119 set(off_side,off_side,get(off_side,off_side)+y);
00120 Complex tmp = get(i,j)-a*y;
00121 set(i,j,tmp);
00122 set(j,i,tmp);
00123 }
00124
00125
00126 void AdmittanceNetwork_PAR1::remove_line(int i, int j, Complex y, double turns,
00127 double phase_shift, int side)
00128 {
00129 Complex a = polar(turns,phase_shift);
00130 set(side,side,get(side,side)-a*a*y);
00131 int off_side = i;
00132 if (side == off_side) off_side = j;
00133 set(off_side,off_side,get(off_side,off_side)-y);
00134 Complex tmp = get(i,j)+a*y;
00135 set(i,j,tmp);
00136 set(j,i,tmp);
00137 }
00138
00139
00140 void AdmittanceNetwork_PAR1::add_self(int i, Complex y)
00141 {
00142 set(i,i,get(i,i)+y);
00143 }
00144
00145
00146 void AdmittanceNetwork_PAR1::solve_for_current(const Complex* V, Complex* I)
00147 {
00148 for (int i = 0; i < data->N; i++)
00149 {
00150 I[i] = Complex(0.0,0.0);
00151 }
00152
00153
00154 SparseComplexMatrix *mat = &data->Y;
00155 SparseComplexMatrix::const_iterator mit;
00156 for (int i = 0; i < data->N; i++)
00157 {
00158 SparseComplexMatrix::const_iterator mit = mat->find( i );
00159 assert( mit != mat->end() );
00160
00161 const SparseComplexVector &vec = mit->second;
00162 for( SparseComplexVector::const_iterator vit = vec.begin();
00163 vit != vec.end(); ++vit )
00164 {
00165 int j = vit->first;
00166 Complex val = vit->second;
00167 I[i] += val*V[j];
00168 }
00169 }
00170 }
00171
00172
00173 static const int
00174 DO_NONE=0, DO_INIT=1, DO_FACTOR=2, DO_SOLVE=3, DO_FINI=4, DO_DONE=5;
00175 #define CMDSTR(c) \
00176 ((c)==DO_NONE?"NONE": \
00177 (c)==DO_INIT?"INIT": \
00178 (c)==DO_FACTOR?"FACTOR": \
00179 (c)==DO_SOLVE?"SOLVE": \
00180 (c)==DO_FINI?"FINI": \
00181 (c)==DO_DONE?"DONE": \
00182 "UNKNOWN")
00183
00184
00185 #define ICNTL(_i) icntl[(_i)-1]
00186
00187
00188 void AdmittanceNetwork_PAR1::par1_init( PAR1Data *par1 )
00189 {
00190 ZMUMPS_STRUC_C *id = &par1->id;
00191
00192
00193 id->par=1; id->sym=0; id->comm_fortran = -987654;
00194
00195
00196 id->irn = 0;
00197 id->jcn = 0;
00198 id->a = 0;
00199 id->rhs = 0;
00200 id->n = 0;
00201
00202
00203 id->job=-1;
00204 zmumps_c(id);
00205
00206 par1->inited = true;
00207 }
00208
00209
00210 void AdmittanceNetwork_PAR1::print_timing( PAR1Data *par1 )
00211 {
00212 cout << "PAR1 TIMINGS " << endl;
00213 cout << "AVERAGES:" <<
00214 " ANALYSIS " << (par1->analysistime/par1->nfactors) <<
00215 " FACTOR " << (par1->factortime/par1->nfactors) <<
00216 " SOLVE " << (par1->solvetime/par1->nsolves) <<
00217 endl;
00218 cout << "TOTALS:" <<
00219 " ANALYSIS * " <<par1->nfactors <<" = " <<par1->analysistime<<
00220 " FACTOR * " <<par1->nfactors << " = " <<par1->factortime<<
00221 " SOLVE * " <<par1->nsolves << " = " <<par1->solvetime<<
00222 endl;
00223 }
00224
00225
00226 void AdmittanceNetwork_PAR1::par1_factor( PAR1Data *par1 )
00227 {
00228 ZMUMPS_STRUC_C *id = &par1->id;
00229
00230
00231 if( false )
00232 {
00233 id->ICNTL(1)=6; id->ICNTL(2)=3; id->ICNTL(3)=6; id->ICNTL(4)=3;
00234 }
00235 else
00236 {
00237 id->ICNTL(1)=0; id->ICNTL(2)=0; id->ICNTL(3)=0; id->ICNTL(4)=0;
00238 }
00239 id->ICNTL(24)=1;
00240 id->ICNTL(7)=6;
00241
00242 par1->nfactors++;
00243 cout << "PAR1 FACTOR STARTED " << par1->nfactors << endl;
00244
00245 double t1 = MPI_Wtime();
00246
00247
00248 id->job=1;
00249 zmumps_c(id);
00250
00251 double t2 = MPI_Wtime();
00252 par1->analysistime += (t2-t1);
00253
00254
00255 id->job=2;
00256 zmumps_c(id);
00257
00258 double t3 = MPI_Wtime();
00259 par1->factortime += (t3-t2);
00260
00261 cout << "PAR1 FACTOR ENDED " << par1->nfactors <<
00262 " INFO " << id->info[0] << endl;
00263 print_timing( par1 );
00264 }
00265
00266
00267 void AdmittanceNetwork_PAR1::par1_solve( PAR1Data *par1 )
00268 {
00269 ZMUMPS_STRUC_C *id = &par1->id;
00270
00271 par1->nsolves++;
00272 cout << "PAR1 SOLVE STARTED " << par1->nsolves << endl;
00273
00274 double t1 = MPI_Wtime();
00275
00276
00277 id->job=3;
00278 zmumps_c(id);
00279
00280 double t2 = MPI_Wtime();
00281 par1->solvetime += (t2-t1);
00282
00283 cout << "PAR1 SOLVE ENDED " << par1->nsolves <<
00284 " INFO " << id->info[0] << endl;
00285 print_timing( par1 );
00286 }
00287
00288
00289 void AdmittanceNetwork_PAR1::par1_fini( PAR1Data *par1 )
00290 {
00291 ZMUMPS_STRUC_C *id = &par1->id;
00292
00293
00294 id->job=-2;
00295 zmumps_c(id);
00296 id->n = 0;
00297
00298 par1->inited = false;
00299 }
00300
00301
00302 void AdmittanceNetwork_PAR1::master_order( PAR1Data *par1, int command )
00303 {
00304 int root = 0;
00305 int err = MPI_Bcast(&command, 1, MPI_INT, root, MPI_COMM_WORLD);
00306 assert( err == 0 );
00307 if( par1 )
00308 {
00309 par1->ncmds++;
00310 cout << "Master cmd " << par1->ncmds << " " << CMDSTR(command) << endl;
00311 }
00312 }
00313
00314
00315 void AdmittanceNetwork_PAR1::master_init( PAR1Data *par1 )
00316 {
00317 master_order( par1, DO_INIT );
00318 par1_init( par1 );
00319 }
00320
00321
00322 static bool doperm = false;
00323 void AdmittanceNetwork_PAR1::master_permute( PAR1Data *par1 )
00324 {
00325 SparseComplexMatrix *mat = &data->Y;
00326
00327 if( par1->nfactors <= 0 )
00328 {
00329 const char *pstr = getenv("DOPERM");
00330 doperm = pstr && !strcmp(pstr,"TRUE");
00331 cout << "DOPERM= " << doperm << endl;
00332 }
00333
00334 if( !doperm )
00335 {
00336 cout<<"+-+-+-+-+-+-+ Skipping Permutation +-+-+-+-+-+-+"<<endl;
00337 return;
00338 }
00339 cout<<"+-+-+-+-+-+-+ Performing Permutation +-+-+-+-+-+-+"<<endl;
00340
00341 BWGraph bwg;
00342 for( int i = 0; i < data->N; i++ ) { bwg.addnode( i ); }
00343
00344 for( SparseComplexMatrix::const_iterator mit = mat->begin();
00345 mit != mat->end(); ++mit )
00346 {
00347 int i = mit->first;
00348 const SparseComplexVector &vec = mit->second;
00349 for( SparseComplexVector::const_iterator citj=vec.begin();
00350 citj != vec.end(); ++citj )
00351 {
00352 int j = citj->first;
00353 bwg.addedge( i, j, false );
00354 }
00355 }
00356
00357 bool dbgprint = true;
00358 if( dbgprint )
00359 {
00360 cout << "#nodes= " << bwg.getnumnodes() << endl;
00361 cout << "#edges= " << bwg.getnumedges() << endl;
00362 Labels identityperm;
00363 initlabels( identityperm, bwg.getnumnodes(), true );
00364 int idbw = bwg.bandwidth(identityperm);
00365 cout << "InputBW= "<< idbw << endl;
00366 }
00367
00368 Labels &perm = data->perm;
00369 perm.clear();
00370 bwg.relabel( perm );
00371
00372 if( dbgprint )
00373 {
00374 int minbw = bwg.bandwidth(perm);
00375 cout << "FinalBW= "<< minbw << endl;
00376 }
00377 }
00378
00379
00380 void AdmittanceNetwork_PAR1::master_factor( PAR1Data *par1 )
00381 {
00382 ZMUMPS_STRUC_C *id = &par1->id;
00383 SparseComplexMatrix *mat = &data->Y;
00384
00385
00386 id->n= data->N;
00387 id->nz= 0;
00388
00389
00390 for( SparseComplexMatrix::const_iterator mit = mat->begin();
00391 mit != mat->end(); ++mit )
00392 {
00393 id->nz += mit->second.size();
00394 }
00395
00396
00397 id->irn = new int[id->nz];
00398 id->jcn = new int[id->nz];
00399 id->a = new ZMUMPS_COMPLEX[id->nz];
00400 id->rhs = new ZMUMPS_COMPLEX[id->n];
00401
00402
00403 master_permute( par1 );
00404
00405
00406 int k = 0;
00407 for( SparseComplexMatrix::const_iterator mit = mat->begin();
00408 mit != mat->end(); ++mit )
00409 {
00410 int i = mit->first;
00411 const SparseComplexVector &vec = mit->second;
00412 for( SparseComplexVector::const_iterator citj=vec.begin();
00413 citj != vec.end(); ++citj )
00414 {
00415 int j = citj->first;
00416 const Complex &val = citj->second;
00417 id->irn[k] = data->perm[i]+1;
00418 id->jcn[k] = data->perm[j]+1;
00419 id->a[k].r = real(val);
00420 id->a[k].i = imag(val);
00421 k++;
00422 }
00423 }
00424 assert( k == id->nz );
00425
00426 master_order( par1, DO_FACTOR );
00427 par1_factor( par1 );
00428 }
00429
00430
00431 void AdmittanceNetwork_PAR1::master_solve( PAR1Data *par1 )
00432 {
00433 master_order( par1, DO_SOLVE );
00434 par1_solve( par1 );
00435 }
00436
00437
00438 void AdmittanceNetwork_PAR1::master_fini( PAR1Data *par1 )
00439 {
00440 ZMUMPS_STRUC_C *id = &par1->id;
00441
00442 master_order( par1, DO_FINI );
00443 par1_fini( par1 );
00444
00445
00446 delete [] id->irn; id->irn = 0;
00447 delete [] id->jcn; id->jcn = 0;
00448 delete [] id->a; id->a = 0;
00449 delete [] id->rhs; id->rhs = 0;
00450 }
00451
00452
00453 static bool verifysol = false;
00454 static bool printsol = false;
00455 void AdmittanceNetwork_PAR1::solve_for_voltage(const Complex* I, Complex* V)
00456 {
00457 ZMUMPS_STRUC_C *id = &data->par1.id;
00458
00459 if( dirty )
00460 {
00461 if( data->par1.inited )
00462 {
00463 master_fini( &data->par1 );
00464 }
00465 master_init( &data->par1 );
00466 master_factor( &data->par1 );
00467 dirty = false;
00468 }
00469
00470
00471 for (int i = 0; i < data->N; i++)
00472 {
00473 int permi = data->perm[i];
00474 id->rhs[permi].r = real(I[i]);
00475 id->rhs[permi].i = imag(I[i]);
00476 }
00477
00478 master_solve( &data->par1 );
00479
00480
00481 for (int i = 0; i < data->N; i++)
00482 {
00483 int permi = data->perm[i];
00484 V[i] = Complex( id->rhs[permi].r, id->rhs[permi].i );
00485 }
00486
00487 if( data->par1.nfactors <= 1 && data->mpi.rank == 0 )
00488 {
00489 const char *pstr = getenv("VERIFYSOLUTIONS");
00490 verifysol = pstr && !strcmp(pstr,"TRUE");
00491 cout << "VERIFYSOLUTIONS= " << verifysol << endl;
00492 }
00493
00494 if( !verifysol )
00495 {
00496 cout<<"+-+-+-+-+-+-+ Skipping verification of solution +-+-+-+-+-+-+"<<endl;
00497 }
00498 else
00499 {
00500 cout<<"+-+-+-+-+-+-+ Performing verification of solution +-+-+-+-+-+-+"<<endl;
00501 double sumreal = 0, sumimag = 0;
00502 Complex *tempI = new Complex[data->N];
00503 solve_for_current( V, tempI );
00504 if(printsol){cout<<"@ "<<data->N<<" [I,tempI]= {";}
00505 for( int i = 0; i < data->N; i++ )
00506 {
00507 #define SQR(x) ((x)*(x))
00508 #define PR(c) ((c)<1e-12?0:(c))
00509 sumreal += SQR(real(I[i])-real(tempI[i]));
00510 sumimag += SQR(imag(I[i])-imag(tempI[i]));
00511 if(printsol){cout<<"\t"<<i<<"\t[ ("<<PR(real(I[i]))<<"=="<<PR(real(tempI[i])) << ") ("<<PR(imag(I[i]))<<"=="<<PR(imag(tempI[i]))<<")"<<endl;}
00512 }
00513 if(printsol){cout<<"}"<<endl;}
00514 double meanreal = sumreal/double(data->N), rmsreal = sqrt(meanreal);
00515 double meanimag = sumimag/double(data->N), rmsimag = sqrt(meanimag);
00516 if( rmsreal > 1e-7 || rmsimag >1e-7 )
00517 {
00518 cerr<<"RMS too high: rmsreal= "<<rmsreal<<" rmsimag= "<<rmsimag << endl;
00519 exit(1);
00520 }
00521 if(printsol){cout << "RMS = " << sqrt(meanreal)<<", "<<sqrt(meanimag) << endl;}
00522 delete [] tempI;
00523 }
00524 }
00525
00526
00527 void AdmittanceNetwork_PAR1::main_begin(int *pac, char ***pav)
00528 {
00529 int rank = 0, nranks = 0;
00530 MPI_Init(pac,pav);
00531 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00532 MPI_Comm_size( MPI_COMM_WORLD, &nranks );
00533
00534 {
00535 char temp[1000]; sprintf(temp,"stdout-%d-%d.txt",nranks,rank);
00536 freopen( temp, "w", stdout );
00537 }
00538
00539 if( rank == 0 )
00540 {
00541
00542 }
00543 else
00544 {
00545 PAR1Data *mdata = new PAR1Data();
00546 while(1)
00547 {
00548 int root = 0, command = DO_NONE;
00549 int err =
00550 MPI_Bcast(&command, 1, MPI_INT, root, MPI_COMM_WORLD);
00551 if( err != 0 || command == DO_DONE )
00552 {
00553 cout << "Slave done " << err << endl;
00554 break;
00555 }
00556 else
00557 {
00558 mdata->ncmds++;
00559 cout << "Slave cmd " << mdata->ncmds << " " << CMDSTR(command) << endl;
00560 switch( command )
00561 {
00562 case DO_INIT:
00563 {
00564 par1_init( mdata );
00565 break;
00566 }
00567 case DO_FACTOR:
00568 {
00569 par1_factor( mdata );
00570 break;
00571 }
00572 case DO_SOLVE:
00573 {
00574 par1_solve( mdata );
00575 break;
00576 }
00577 case DO_FINI:
00578 {
00579 par1_fini( mdata );
00580 break;
00581 }
00582 default:
00583 {
00584 cout << "Slave unknown cmd " << command << endl;
00585 assert( false );
00586 break;
00587 }
00588 }
00589 }
00590 }
00591 delete mdata;
00592 main_end( rank );
00593 if( rank != 0 )
00594 {
00595 cout << "Slave exit" << endl;
00596 exit(0);
00597 }
00598 }
00599 }
00600
00601
00602 void AdmittanceNetwork_PAR1::main_end( int rank )
00603 {
00604 cout << "main_end begin" << endl;
00605 if( rank == 0 )
00606 {
00607 master_order( 0, DO_DONE );
00608 }
00609 cout << "MPI_Barrier" << endl;
00610 MPI_Barrier(MPI_COMM_WORLD);
00611 cout << "MPI_Finalize" << endl;
00612 MPI_Finalize();
00613 cout << "main_end end" << endl;
00614 }
00615
00616