Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_MC64.cpp
Go to the documentation of this file.
1#include "Amesos_ConfigDefs.h"
2#include "Amesos_MC64.h"
3#include "Epetra_Comm.h"
4#include "Epetra_RowMatrix.h"
5#include <map>
6
7extern "C" void F77_FUNC(mc64id, MC64ID)(int*);
8extern "C" void F77_FUNC(mc64ad, MC64AD)(int*, int*, int*, int*, int*,
9 double*, int*, int*, int*, int*,
10 int*, double*, int*, int*);
11
12// ===========================================================================
13Amesos_MC64::Amesos_MC64(const Epetra_RowMatrix& A, int JOB,
14 const bool StoreTranspose, const bool analyze) :
15 A_(A)
16{
17 if (A_.Comm().NumProc() != 1)
18 {
19 std::cerr << "Class Amesos_MC64 can be used with one processor only!" << std::endl;
20 exit(EXIT_FAILURE);
21 }
22 F77_FUNC(mc64id, MC64ID)(ICNTL_);
23
24 Compute(JOB, StoreTranspose, analyze);
25}
26
27// ===========================================================================
28int Amesos_MC64::Compute(int JOB, const bool StoreTranspose, const bool analyze)
29{
30 // convert A_ into column-based format
31
32 int MaxNumEntries = A_.MaxNumEntries();
33 int N = A_.NumMyRows();
34 int NE = A_.NumMyNonzeros();
35 std::vector<int> IP;
36 std::vector<int> IRN;
37 std::vector<double> VAL;
38
39 std::vector<int> Indices(MaxNumEntries);
40 std::vector<double> Values(MaxNumEntries);
41
42 if (StoreTranspose)
43 {
44 // cheapest way, just store the transpose of A and not A. This is because
45 // we have easy access to rows and not to columns.
46 IP.resize(N + 1); IP[0] = 1;
47 IRN.resize(NE);
48 VAL.resize(NE);
49 int count = 0;
50
51 for (int i = 0 ; i < N ; ++i)
52 {
53 int NumEntries = 0;
54
55 A_.ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0],
56 &Indices[0]);
57
58 IP[i + 1] = IP[i] + NumEntries;
59
60 for (int j = 0 ; j < NumEntries ; ++j)
61 {
62 IRN[count] = Indices[j] + 1;
63 VAL[count] = Values[j];
64 ++count;
65 }
66 }
67 assert(count == NE);
68 }
69 else
70 {
71 // this stores the matrix and not its transpose, but it is more memory
72 // demading. The ifdef'd out part is simple and fast, but very memory
73 // demanding.
74
75#if 0
76 IRN.resize(N * MaxNumEntries);
77 VAL.resize(N * MaxNumEntries);
78
79 std::vector<int> count(N);
80 for (int i = 0 ; i < N ; ++i) count[i] = 0;
81
82 for (int i = 0 ; i < N ; ++i)
83 {
84 int NumEntries = 0;
85
86 A_.ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0],
87 &Indices[0]);
88
89 for (int j = 0 ; j < NumEntries ; ++j)
90 {
91 int col = Indices[j];
92 IRN[col * MaxNumEntries + count[col]] = i + 1;
93 VAL[col * MaxNumEntries + count[col]] = Values[j];
94 ++count[col];
95 }
96 }
97
98 // now compact storage
99 int k = 0;
100 for (int col = 0 ; col < N ; ++col)
101 {
102 for (int row = 0 ; row < count[col] ; ++row)
103 {
104 IRN[k] = IRN[col * MaxNumEntries + row];
105 VAL[k] = VAL[col * MaxNumEntries + row];
106 ++k;
107 }
108 }
109 assert (k == NE);
110
111 IRN.resize(k);
112 VAL.resize(k);
113
114 IP.resize(N + 1);
115 IP[0] = 1;
116
117 for (int col = 0 ; col < N ; ++col)
118 IP[col + 1] = IP[col] + count[col];
119#else
120 std::vector<std::vector<int> > cols(N);
121 std::vector<std::vector<double> > vals(N);
122
123 for (int i = 1 ; i <= N ; ++i)
124 {
125 int NumEntries = 0;
126
127 A_.ExtractMyRowCopy(i - 1, MaxNumEntries, NumEntries, &Values[0],
128 &Indices[0]);
129
130 for (int j = 0 ; j < NumEntries ; ++j)
131 {
132 cols[Indices[j]].push_back(i);
133 vals[Indices[j]].push_back(Values[j]);
134 }
135 }
136
137 IP.resize(N + 1); IP[0] = 1;
138 IRN.resize(NE);
139 VAL.resize(NE);
140 int count = 0;
141
142 for (int i = 0 ; i < N ; ++i)
143 {
144 IP[i + 1] = IP[i] + cols[i].size();
145
146 for (int j = 0 ; j < cols[i].size() ; ++j)
147 {
148 IRN[count] = cols[i][j];
149 VAL[count] = vals[i][j];
150 ++count;
151 }
152 }
153#endif
154 }
155
156 int NUM;
157 CPERM_.resize(N);
158 int LIW = 10 * N + NE;
159 std::vector<int> IW(LIW);
160 int LDW = 3 * N + NE;
161 DW_.resize(LDW);
162
163 JOB = 5;
164 F77_FUNC(mc64ad, MC64aD)(&JOB, &N, &NE, &IP[0], &IRN[0], &VAL[0],
165 &NUM, &CPERM_[0], &LIW, &IW[0], &LDW,
166 &DW_[0], ICNTL_, INFO_);
167
168 if (analyze)
169 {
170 std::map<double, int> table;
171 for (int col = 0 ; col < N ; ++col)
172 {
173 for (int j = IP[col] ; j < IP[col + 1] ; ++j)
174 {
175 int row = IRN[j - 1] - 1;
176 int new_col = CPERM_[col];
177 double new_val = VAL[j - 1] * exp(DW_[row] + DW_[col + N]);
178 if (new_val < 0.0) new_val = -new_val;
179 if (new_val > 0.1) table[0.1]++;
180 else if (new_val > 0.01) table[0.01]++;
181 else if (new_val > 0.001) table[0.001]++;
182 else if (new_val > 0.0001) table[0.0001]++;
183 else if (new_val > 0.00001) table[0.00001]++;
184 else if (new_val > 0.000001) table[0.000001]++;
185 else table[0.0]++;
186 }
187 }
188
189 std::cout << "# elements (total) = " << NE << std::endl;
190 std::cout << "# elements > 0.1 = " << table[0.1] << std::endl;
191 std::cout << "# elements > 0.01 = " << table[0.01] << std::endl;
192 std::cout << "# elements > 0.001 = " << table[0.001] << std::endl;
193 std::cout << "# elements > 0.0001 = " << table[0.0001] << std::endl;
194 std::cout << "# elements > 0.00001 = " << table[0.00001] << std::endl;
195 std::cout << "# elements > 0.000001 = " << table[0.000001] << std::endl;
196 std::cout << "# elements <=0.000001 = " << table[0.0] << std::endl;
197 }
198
199 AMESOS_RETURN(INFO_[0]);
200}
201
202// ===========================================================================
203double* Amesos_MC64::GetColScaling()
204{
205 return((double*)&DW_[0 + A_.NumMyRows()]);
206}
#define AMESOS_RETURN(amesos_err)
void F77_FUNC(mc64id, MC64ID)(int *)