#include "ml_include.h"
#if defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS) && defined(HAVE_ML_EPETRAEXT) && defined(HAVE_ML_AZTECOO)
#ifdef HAVE_MPI
#include "mpi.h"
#endif
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
#include "Epetra_Export.h"
#include "AztecOO.h"
#include "ml_epetra.h"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_DefaultComm.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "EpetraExt_RowMatrixOut.h"
#include "EpetraExt_VectorOut.h"
#include "EpetraExt_BlockMapOut.h"
#include "EpetraExt_BlockMapIn.h"
#include "EpetraExt_CrsMatrixIn.h"
#include "EpetraExt_VectorIn.h"
int MatrixMarketFileToCrsMatrix(const char *filename,
{
return(EpetraExt::MatrixMarketFileToCrsMatrixHandle(filename, rowMap.
Comm(), A,
&rowMap,NULL,
&rangeMap, &domainMap));
}
return EpetraExt::MatrixMarketFileToVector(filename,map,v);
}
ML_Comm *mlcomm;
int main(int argc, char *argv[])
{
#ifdef ML_MPI
MPI_Init(&argc,&argv);
int commWorldSize;
MPI_Comm_size(MPI_COMM_WORLD,&commWorldSize);
std::vector<int> splitKey;
int rankToDrop = 1;
for (int i=0; i<commWorldSize; ++i)
splitKey.push_back(0);
splitKey[rankToDrop] = MPI_UNDEFINED;
int myrank;
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm subcomm;
MPI_Comm_split(MPI_COMM_WORLD, splitKey[myrank], myrank, &subcomm);
if (myrank == rankToDrop) goto droppedRankLabel;
#endif
{
#ifdef ML_MPI
#else
#endif
ML_Comm_Create(&mlcomm);
char *datafile;
#ifdef CurlCurlAndMassAreSeparate
if (argc < 5 && argc > 9) {
std::cout << "usage: ml_maxwell.exe <S> <M> <T> <Kn> [rhs] [xml file] [edge map] [node map]"
<< std::endl;
std::cout << " S = edge stiffness matrix file" << std::endl;
std::cout << " M = edge mass matrix file" << std::endl;
std::cout << " T = discrete gradient file" << std::endl;
std::cout << " Kn = auxiliary nodal FE matrix file" << std::endl;
std::cout << " rhs = rhs vector" << std::endl;
std::cout << " xml file = xml solver options" <<std::endl;
std::cout << " edge map = edge distribution over processors" << std::endl;
std::cout << " node map = node distribution over processors" << std::endl;
std::cout << argc << std::endl;
}
#else
if (argc < 4 && argc > 8) {
std::cout << "usage: ml_maxwell.exe <A> <T> <Kn> [rhs] [xml file] [edge map] [node map]" <<std::endl;
std::cout << " A = edge element matrix file" << std::endl;
std::cout << " T = discrete gradient file" << std::endl;
std::cout << " Kn = auxiliary nodal FE matrix file" << std::endl;
std::cout << " rhs = rhs vector" << std::endl;
std::cout << " xml file = xml solver options" <<std::endl;
std::cout << " edge map = edge distribution over processors" << std::endl;
std::cout << " node map = node distribution over processors" << std::endl;
std::cout << argc << std::endl;
}
#endif
#ifdef ML_MPI
MPI_Finalize();
#endif
exit(1);
}
#ifdef CurlCurlAndMassAreSeparate
if (argc > 7)
#else
if (argc > 6)
#endif
{
datafile = argv[7];
printf("Reading in edge map from %s ...\n",datafile);
fflush(stdout);
}
EpetraExt::MatrixMarketFileToMap(datafile, Comm, edgeMap);
datafile = argv[6];
printf("Reading in node map from %s ...\n",datafile);
fflush(stdout);
}
EpetraExt::MatrixMarketFileToMap(datafile, Comm, nodeMap);
}
else {
printf("Using linear edge and node maps ...\n");
const int lineLength = 1025;
char line[lineLength];
FILE *handle;
int M,N,NZ;
#ifdef CurlCurlAndMassAreSeparate
handle = fopen(argv[3],"r");
#else
handle = fopen(argv[2],"r");
#endif
if (handle == 0) EPETRA_CHK_ERR(-1);
do {
if(fgets(line, lineLength, handle)==0) {if (handle!=0) fclose(handle);}
} while (line[0] == '%');
if(sscanf(line,"%d %d %d", &M, &N, &NZ)==0) {if (handle!=0) fclose(handle);}
fclose(handle);
}
Teuchos::ParameterList MLList;
int *options = new int[AZ_OPTIONS_SIZE];
double *params = new double[AZ_PARAMS_SIZE];
#ifdef CurlCurlAndMassAreSeparate
const int xml_idx = 6;
#else
const int xml_idx = 5;
#endif
if (argc > xml_idx) {
Teuchos::updateParametersFromXmlFileAndBroadcast(argv[xml_idx],Teuchos::Ptr<Teuchos::ParameterList>(&MLList),*Teuchos::DefaultComm<int>::getComm());
}
else {
MLList.set("ML output", 10);
MLList.set("aggregation: type", "Uncoupled");
MLList.set("coarse: max size", 128);
MLList.set("aggregation: threshold", 0.0);
MLList.set("subsmoother: type", "symmetric Gauss-Seidel");
MLList.set("max levels", 2);
MLList.set("aggregation: damping factor",0.0);
MLList.set("coarse: type", "Amesos-KLU");
MLList.set("aggregation: do qr", false);
MLList.set("smoother: sweeps",1);
MLList.set("subsmoother: edge sweeps",1);
MLList.set("subsmoother: node sweeps",1);
MLList.set("smoother: Hiptmair efficient symmetric",false);
}
#ifdef CurlCurlAndMassAreSeparate
for (int i = 1; i <6; i++) {
datafile = argv[i];
printf("reading %s ....\n",datafile); fflush(stdout);
}
switch (i) {
case 1:
MatrixMarketFileToCrsMatrix(datafile,*edgeMap,*edgeMap,*edgeMap,CurlCurl);
break;
case 2:
MatrixMarketFileToCrsMatrix(datafile, *edgeMap, *edgeMap, *edgeMap, Mass);
break;
case 3:
MatrixMarketFileToCrsMatrix(datafile, *edgeMap, *edgeMap,*nodeMap, T);
break;
case 4:
MatrixMarketFileToCrsMatrix(datafile, *nodeMap,*nodeMap, *nodeMap, Kn);
break;
case 5:
MatrixMarketFileToVector(datafile, *edgeMap,rhs);
break;
}
}
#else
for (int i = 1; i <5; i++) {
datafile = argv[i];
printf("reading %s ....\n",datafile); fflush(stdout);
}
switch (i) {
case 1:
MatrixMarketFileToCrsMatrix(datafile,*edgeMap,*edgeMap,*edgeMap,CCplusM);
break;
case 2:
MatrixMarketFileToCrsMatrix(datafile, *edgeMap, *edgeMap, *nodeMap, T);
break;
case 3:
MatrixMarketFileToCrsMatrix(datafile, *nodeMap, *nodeMap, *nodeMap, Kn);
break;
case 4:
MatrixMarketFileToVector(datafile, *edgeMap,rhs);
break;
}
}
#endif
std::cout<<"*** ML Parameters ***\n"<<MLList<<std::endl;
#ifdef CurlCurlAndMassAreSeparate
#endif
#ifdef CurlCurlAndMassAreSeparate
#else
#endif
if(!rhs) {
std::cout << "Putting in a zero initial guess and random rhs (in the range of S+M)" << std::endl;
x.Random();
}
else {
std::cout << "Using user rhs"<<std::endl;
}
x.PutScalar(0.0);
double vecnorm;
if (Comm.
MyPID() == 0) std::cout <<
"||rhs|| = " << vecnorm << std::endl;
x.Norm2(&vecnorm);
if (Comm.
MyPID() == 0) std::cout <<
"||x|| = " << vecnorm << std::endl;
AztecOO solver(Problem);
solver.SetPrecOperator(MLPrec);
solver.SetAztecOption(AZ_solver, AZ_cg);
solver.SetAztecOption(AZ_output, 1);
solver.Iterate(15, 1e-10);
delete MLPrec;
delete CurlCurl;
delete CCplusM;
delete Mass;
delete T;
delete Kn;
delete rhs;
delete nodeMap;
delete edgeMap;
delete [] params;
delete [] options;
ML_Comm_Destroy(&mlcomm);
}
droppedRankLabel:
#ifdef ML_MPI
MPI_Finalize();
#endif
return 0;
}
#else
#include <stdlib.h>
#include <stdio.h>
#ifdef HAVE_MPI
#include "mpi.h"
#endif
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
#endif
puts("Please configure ML with:");
#if !defined(HAVE_ML_EPETRA)
puts("--enable-epetra");
#endif
#if !defined(HAVE_ML_TEUCHOS)
puts("--enable-teuchos");
#endif
#if !defined(HAVE_ML_EPETRAEXT)
puts("--enable-epetraext");
#endif
#if !defined(HAVE_ML_AZTECOO)
puts("--enable-aztecoo");
#endif
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return 0;
}
#endif
const Epetra_Comm & Comm() const
const Epetra_Map & DomainMap() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
int Norm2(double *Result) const
MultiLevelPreconditioner: a class to define black-box multilevel preconditioners using aggregation me...
Definition ml_MultiLevelPreconditioner.h:241
void PrintUnused() const
Prints unused parameters in the input ParameterList on standard output.
Definition ml_MultiLevelPreconditioner.h:406
void Print(int level=-2)
Print the individual operators in the multigrid hierarchy.
Definition ml_MultiLevelPreconditioner.cpp:2807
Interface to the Trilinos package Anasazi.
Epetra_CrsMatrix * Epetra_MatrixAdd(Epetra_RowMatrix *B, Epetra_RowMatrix *Bt, double scalar)
Adds two Epetra_RowMatrix's, returns the result as an Epetra_CrsMatrix.
int SetDefaults(std::string ProblemType, Teuchos::ParameterList &List, int *options=0, double *params=0, const bool OverWrite=true)
Sets default parameters for aggregation-based preconditioners.