examples/ml_maxwell_test.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002     ml_maxwell_test.cpp  -  test program for Trilinos ML in Maxwell mode
00003                              -------------------
00004     begin                : Tue Apr 27 2004
00005     copyright            : (C) 2004 by Roman Geus
00006     email                : roman.geus@psi.ch
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 // #include "Epetra_config.h"
00019 #include <ml_config.h>
00020 
00021 #ifdef HAVE_MPI
00022 #include <mpi.h>
00023 #include <Epetra_MpiComm.h>
00024 #else
00025 #include <Epetra_SerialComm.h>
00026 #endif
00027 
00028 // for the ML Maxwell preconditioner...
00029 #if !defined(HAVE_ML_EPETRA) || !defined(HAVE_ML_TEUCHOS)
00030 #error Trilinos needs to be configured with --enable-epetra --enable-teuchos
00031 #endif
00032 
00033 #include <ml_include.h>
00034 #include <ml_epetra_preconditioner.h>
00035 #include <Teuchos_ParameterList.hpp>
00036 
00037 #include <AztecOO.h>
00038 #include <Epetra_Operator.h>
00039 #include <Epetra_MultiVector.h>
00040 #include <Epetra_Vector.h>
00041 #include <Epetra_LinearProblem.h>
00042 #include "Epetra_Export.h"
00043 
00044 #include <iostream>
00045 #include <fstream>
00046 #include <cstdlib>
00047 #include <unistd.h>
00048 
00049 #include "matrixmarket/mtxreader0.h"
00050 #include "matrixmarket/mtxreader.h"
00051 
00052 void finalize() {
00053 #ifdef HAVE_MPI
00054     MPI_Finalize();
00055 #endif
00056 }
00057 
00058 int permute(const Epetra_Map& newHMap,
00059             const Epetra_Map& newKMap,
00060             Epetra_CrsMatrix*& K, 
00061             Epetra_CrsMatrix*& H, 
00062             Epetra_CrsMatrix*& Y) {
00063 
00064     Epetra_Export HExport(H->RangeMap(), newHMap);
00065     Epetra_CrsMatrix *Hnew = new Epetra_CrsMatrix(Copy, newHMap, 0);
00066     Hnew->Export(*H, HExport, Insert);
00067     Hnew->FillComplete(newHMap, newHMap);
00068     delete H;
00069     H = Hnew;
00070 
00071     Epetra_Export KExport(K->RangeMap(), newKMap);
00072     Epetra_CrsMatrix *Knew = new Epetra_CrsMatrix(Copy, newKMap, 0);
00073     Knew->Export(*K, KExport, Insert);
00074     Knew->FillComplete(newKMap, newKMap);
00075     delete K;
00076     K = Knew;
00077 
00078     Epetra_Export YExport(Y->RangeMap(), newKMap);
00079     Epetra_CrsMatrix *Ynew = new Epetra_CrsMatrix(Copy, newKMap, 0);
00080     Ynew->Export(*Y, YExport, Insert);
00081     Ynew->FillComplete(H->DomainMap(), newKMap);
00082     delete Y;
00083     Y = Ynew;
00084 
00085     return(0);
00086 }
00087 
00088 int permute_ml(Epetra_CrsMatrix*& K, Epetra_CrsMatrix*& H, Epetra_CrsMatrix*& Y) {
00089     const Epetra_Comm& Comm(K->Comm());
00090     const Epetra_Map& newHMap = H->RangeMap();
00091     const Epetra_Map& newMMap = K->RangeMap();
00092 
00093     int num_my_edges = newMMap.NumMyElements();
00094     int num_my_nodes = newHMap.NumMyElements();
00095     
00096     Epetra_Export YExport(Y->RangeMap(), newMMap);
00097     Epetra_CrsMatrix *Ynew = new Epetra_CrsMatrix(Copy, newMMap, 0);
00098     Ynew->Export(*Y, YExport, Insert);
00099     Ynew->FillComplete(newHMap, newMMap);
00100     delete Y;
00101     Y = Ynew;
00102 
00103     int num_my_elements = num_my_edges;
00104     int my_pid = Comm.MyPID();
00105     int num_proc = Comm.NumProc();
00106 
00107     int* local_size = new int[num_proc];
00108     Comm.GatherAll(&num_my_elements, local_size, 1);
00109 
00110     int i = 0;
00111     int j = 0;
00112     int maxs = 0;
00113     for (i = 0; i < num_proc; i ++)
00114         if (local_size[i] > maxs) 
00115             maxs = local_size[i]; 
00116 
00117     int* global_element = new int[maxs * num_proc];
00118     int* next_location = new int[maxs * num_proc];
00119     int offset = my_pid * maxs;
00120 
00121     int* col = new int[2];
00122     int* proc = new int[2];
00123     int* loc = new int[2];
00124     int num_entries = 0;
00125     int *indices = 0;
00126     double *values = 0;
00127     bool notfound = false;
00128 
00129     Epetra_Util util;
00130     util.SetSeed(763206768);
00131 
00132     //
00133     //  my replacement ("location") of the Epetra_Map.RemoteIDList method
00134     //
00135     int* local_s = new int[num_proc];
00136     local_s[my_pid] = num_my_nodes;
00137     Comm.GatherAll(local_s + my_pid, local_s, 1);
00138     int maxhs = 0;
00139     for (i = 0; i < num_proc; i ++)
00140         if (local_s[i] > maxhs) 
00141             maxhs = local_s[i];
00142     int* global_e = new int[maxhs * num_proc];
00143     int* location = new int[newHMap.NumGlobalElements()];
00144     int offs = my_pid * maxhs;
00145     for (i = 0; i < num_my_nodes; i ++) {
00146         global_e[offs + i] = newHMap.GID(i);
00147     }
00148     Comm.GatherAll(global_e + offs, global_e, maxhs);
00149     for (i = 0; i < num_proc; i++)
00150         for (j = i * maxhs; j < i * maxhs + local_s[i]; j ++) 
00151             location[global_e[j]] = i; 
00152     delete[] local_s;
00153     delete[] global_e;
00154     //
00155     // end 
00156     //
00157 
00158     for (i = 0; i < num_my_elements; i ++) {
00159         global_element[offset + i] = newMMap.GID(i);
00160         assert( Y->ExtractMyRowView(i, num_entries, values, indices) >= 0 );
00161         j = 0;
00162         notfound = true;
00163         while ((j < num_entries) && notfound) {
00164             col[j] = Y->GCID(indices[j]);
00165 
00166             //            newHMap.RemoteIDList(1, col + j, proc + j, loc + j);
00167             proc[j] = location[col[j]];
00168 
00169             if (proc[j] == my_pid) 
00170                 notfound = false;
00171             else
00172                 j ++;
00173         }
00174         if (notfound && (num_entries != 0)) {
00175             if (num_entries == 1)
00176                 j = 0;
00177             else 
00178                 j = util.RandomInt() % 2;
00179             next_location[offset + i] = proc[j];
00180         }
00181         else {
00182             next_location[offset + i] = my_pid; 
00183         }
00184     }
00185     for (i = offset + num_my_elements; i < offset + maxs; i ++)
00186         next_location[i] = -1; 
00187  
00188     Comm.GatherAll(next_location + offset, next_location, maxs);
00189     Comm.GatherAll(global_element + offset, global_element, maxs);    
00190 
00191     num_my_elements = 0;
00192     for (i = 0; i < maxs*num_proc; i ++) {
00193         if ( next_location[i] == my_pid )
00194             global_element[num_my_elements++] = global_element[i];
00195     }
00196 
00197     Epetra_Map MMap(newMMap.NumGlobalElements(), num_my_elements, global_element, 0, Comm);
00198 
00199     delete[] local_size;
00200     delete[] global_element;
00201     delete[] next_location;
00202     delete[] col;
00203     delete[] proc;
00204     delete[] loc;
00205     delete[] location;
00206 
00207     permute(newHMap, MMap, K, H, Y);
00208 
00209     return(0);
00210 }
00211 
00212 
00213 int main(int argc, char *argv[])
00214 {
00215 #ifdef HAVE_MPI
00216     MPI_Init(&argc, &argv);
00217     Epetra_MpiComm Comm(MPI_COMM_WORLD);
00218 #else
00219     Epetra_SerialComm Comm;
00220 #endif
00221 
00222 #ifdef __GNUC__
00223     std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
00224 #endif
00225 
00226     if (Comm.MyPID() == 0)
00227         cout << "running on " << Comm.NumProc() << " processor(s)..." << endl;
00228 
00229     // register atexit function
00230     atexit(finalize);
00231 
00232     //
00233     // read and distribute MatrixMarket matrices
00234     //
00235     BaseMtxReader *reader = 0;
00236 
00237     // edge matrix
00238     string filename("K.mtx");
00239     if (access(filename.c_str(), R_OK)) {
00240         cerr << "Unable to access edge matrix file.\n";
00241         exit(1);
00242     }
00243     reader = new MtxReader0(filename, Comm);
00244     Epetra_CrsMatrix* pK = reader->read();
00245     delete reader;
00246 
00247     // node matrix
00248     filename = "H.mtx";
00249     if (access(filename.c_str(), R_OK)) {
00250         cerr << "Unable to access node matrix file.\n";
00251         exit(1);
00252     }
00253     reader = new MtxReader0(filename, Comm);
00254     Epetra_CrsMatrix* pH = reader->read();
00255     delete reader;
00256 
00257     // transfer matrix (with compatible maps)
00258     const Epetra_Map& edge_map = pK->DomainMap();
00259     const Epetra_Map& node_map = pH->RangeMap();
00260     filename = "Y.mtx";
00261     if (access(filename.c_str(), R_OK)) {
00262         cerr << "Unable to access transfer matrix file.\n";
00263         exit(1);
00264     }
00265     reader = new MtxReader(filename, edge_map, node_map, edge_map);
00266     Epetra_CrsMatrix* pY = reader->read();
00267     delete reader;
00268 
00269     // Permute matrices to make sure that edges are assigned to a processor, which also
00270     // holds at least one of its nodes.
00271     permute_ml(pK, pH, pY);
00272 
00273     //
00274     // Construct ML Maxwell preconditioner
00275     //
00276     Epetra_Operator *prec = 0;
00277     Teuchos::ParameterList MLList;
00278 
00279 #ifndef SKIP_ML
00280     // set defaults for maxwell type problems
00281     ML_Epetra::SetDefaults("maxwell", MLList);
00282     MLList.set("aggregation: type", "Uncoupled-MIS");
00283         
00284     // Hack to bypass a typo in ML_Epetra::SetDefaults
00285     MLList.set("coarse: type", "Amesos-KLU");
00286 
00287     if (Comm.MyPID() == 0)
00288         MLList.print(cout, 8);
00289     
00290     // construct ML preconditioner
00291     prec = new ML_Epetra::MultiLevelPreconditioner(*pK,
00292                                                    *pY,
00293                                                    *pH,
00294                                                    MLList,
00295                                                    true);
00296 #endif
00297     //
00298     // Solve linear system K * x = b
00299     //
00300 
00301     // Right hand side
00302     Epetra_Vector b(pK->OperatorRangeMap());
00303     b.PutScalar(1.0);
00304 
00305     // Initial guess
00306     Epetra_Vector x(pK->OperatorDomainMap());
00307     x.PutScalar(0.0);
00308 
00309     // Setup iterative solver
00310     Epetra_LinearProblem problem(pK, &x, &b);
00311     AztecOO solver(problem);
00312     solver.SetAztecOption(AZ_solver, AZ_bicgstab);
00313 
00314 #ifndef SKIP_ML
00315     solver.SetPrecOperator(prec);
00316 #else
00317     solver.SetAztecOption(AZ_precond, AZ_none);
00318 #endif
00319 
00320     // Solve
00321     solver.Iterate(1000, 1.0E-8);
00322 
00323     // Compute Residual
00324     Epetra_Vector bcomp(pK->OperatorRangeMap());
00325     pK->Apply(x, bcomp);
00326     Epetra_Vector resid(pK->OperatorRangeMap());
00327     resid.Update(1.0, b, -1.0, bcomp, 0.0);
00328     double residual;
00329     resid.Norm2(&residual);
00330     if (Comm.MyPID() == 0)
00331         cout << "Residual    = " << residual << endl;
00332 
00333     // Free memory
00334     delete prec;
00335     delete pK;
00336     delete pH;
00337     delete pY;
00338 
00339     return EXIT_SUCCESS;
00340 }

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7