Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
example_MC64.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Interface to Direct Solver Libraries
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29#include "Amesos_ConfigDefs.h"
30
31#ifdef HAVE_MPI
32#include "Epetra_MpiComm.h"
33#else
34#include "Epetra_SerialComm.h"
35#endif
36#include "Epetra_CrsMatrix.h"
37#include "Epetra_Map.h"
38#include "Amesos_MC64.h"
39
41{
42 public:
43 Ifpack_ReorderOperator(const int size, int* perm, double* scale) :
44 size_(size),
45 perm_(perm),
46 scale_(scale)
47 { }
48
49 int Apply(double* x, double* y)
50 {
51 for (int i = 0 ; i < size_ ; ++i)
52 y[i] = scale_[i] * x[perm_[i]];
53 }
54 private:
55 const int size_;
56 int* perm_;
57 double* scale_;
58};
59
60// Creates the following matrix:
61//
62// | 3.00 5.00 |
63// | 2.00 3.00 0.00 |
64// | 3.00 4.00 6.00 |
65// | 1.00 2.00 |
66//
67// as reported in the HSL MC64 version 1.3.0 user's guide. This simple driver
68// calls MC64 using JOB=5, then prints on screen the CPERM and DW vectors. The
69// reordered and scaled matrix looks like:
70//
71// | 1.00 1.00 |
72// | 0.90 1.00 0.00 |
73// | 1.00 1.00 1.00 |
74// | 0.75 1.00 |
75//
76// \warning The interface between Amesos and MC64 works for serial
77// computations only.
78//
79// \author Marzio Sala, ETHZ.
80//
81// \date Last updated on 05-Feb-06.
82
83int main(int argc, char *argv[])
84{
85 // initialize MPI and Epetra communicator
86#ifdef HAVE_MPI
87 MPI_Init(&argc,&argv);
88 Epetra_MpiComm Comm(MPI_COMM_WORLD);
89#else
90 Epetra_SerialComm Comm;
91#endif
92
93 if (Comm.NumProc() != 1)
94 {
95 std::cerr << "This example can be run with one processor only!" << std::endl;
96#ifdef HAVE_MPI
97 MPI_Finalize();
98#endif
99 exit(EXIT_SUCCESS);
100 }
101
102 int NumMyRows = 4;
103 Epetra_Map Map(NumMyRows, 0, Comm);
104 Epetra_CrsMatrix A(Copy, Map, 0);
105
106 int col;
107 double val;
108
109 col = 0; val = 3.0;
110 A.InsertGlobalValues(0, 1, &val, &col);
111
112 col = 1; val = 5.0;
113 A.InsertGlobalValues(0, 1, &val, &col);
114
115 col = 0; val = 2.0;
116 A.InsertGlobalValues(1, 1, &val, &col);
117
118 col = 1; val = 3.0;
119 A.InsertGlobalValues(1, 1, &val, &col);
120
121 col = 3; val = 0.0;
122 A.InsertGlobalValues(1, 1, &val, &col);
123
124 col = 0; val = 3.0;
125 A.InsertGlobalValues(2, 1, &val, &col);
126
127 col = 2; val = 4.0;
128 A.InsertGlobalValues(2, 1, &val, &col);
129
130 col = 3; val = 6.0;
131 A.InsertGlobalValues(2, 1, &val, &col);
132
133 col = 2; val = 1.0;
134 A.InsertGlobalValues(3, 1, &val, &col);
135
136 col = 3; val = 2.0;
137 A.InsertGlobalValues(3, 1, &val, &col);
138
139 A.FillComplete();
140
141 // Creates an instance of the MC64 reordering and scaling interface
142 // and computes the (column) reordering and scaling, using JOB=5
143 Amesos_MC64 MC64(A, 5);
144
145 // checks the return value
146 std::cout << "INFO(1) = " << MC64.GetINFO(1) << std::endl;
147
148 // Gets the pointer to reordering (CPERM) and scaling (DW). Both
149 // vectors are allocated and free'd within Amesos_MC64.
150 int* CPERM = MC64.GetCPERM();
151 double* DW = MC64.GetDW();
152
153 for (int i = 0 ; i < A.NumMyRows() ; ++i)
154 std::cout << "CPERM[" << i << "] = " << CPERM[i] << std::endl;
155
156 for (int i = 0 ; i < A.NumMyRows() * 2 ; ++i)
157 std::cout << "DW[" << i << "] = " << DW[i] << std::endl;
158
159
160 Ifpack_ReorderOperator RowPerm(4, MC64.GetRowPerm(), MC64.GetRowScaling());
161 Ifpack_ReorderOperator ColPerm(4, MC64.GetColPerm(), MC64.GetColScaling());
162
163#ifdef HAVE_MPI
164 MPI_Finalize() ;
165#endif
166
167 return(EXIT_SUCCESS);
168}
Interface to MC64, reordering and scaling algorithm.
Ifpack_ReorderOperator(const int size, int *perm, double *scale)
int Apply(double *x, double *y)
int main(int argc, char *argv[])