Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
CreateTridi.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Direct Sparse Solver Package
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//
30// CreateTridi populates an empty EpetraCrsMatrix with a tridiagonal with
31// -1 on the off-diagonals and 2 on the diagonal.
32//
33// CreateTridiPlus creates the same matrix as CreateTridi except that it adds
34// -1 in the two off diagonal corners.
35//
36// This code was plaguerized from epetra/example/petra_power_method/cxx_main.cpp
37// presumably written by Mike Heroux.
38//
39// Adapted by Ken Stanley - Aug 2003
40//
41#include "Epetra_Map.h"
42#include "Epetra_CrsMatrix.h"
43
44int CreateTridi(Epetra_CrsMatrix& A)
45{
46
47 Epetra_Map Map = A.RowMap();
48 int NumMyElements = Map.NumMyElements();
49 int NumGlobalElements = Map.NumGlobalElements();
50
51 int * MyGlobalElements = new int[NumMyElements];
52 Map.MyGlobalElements(MyGlobalElements);
53
54 // Add rows one-at-a-time
55 // Need some vectors to help
56 // Off diagonal Values will always be -1
57
58
59 double *Values = new double[3];
60 int *Indices = new int[3];
61#ifndef NDEBUG
62 int NumEntries;
63#endif
64 for (int i=0; i<NumMyElements; i++)
65 {
66 if (MyGlobalElements[i]==0)
67 {
68 Indices[0] = 0;
69 Indices[1] = 1;
70 Values[0] = 2.0;
71 Values[1] = -1.0;
72#ifndef NDEBUG
73 NumEntries = 2;
74#endif
75 }
76 else if (MyGlobalElements[i] == NumGlobalElements-1)
77 {
78 Indices[0] = NumGlobalElements-1;
79 Indices[1] = NumGlobalElements-2;
80 Values[0] = 2.0;
81 Values[1] = -1.0;
82#ifndef NDEBUG
83 NumEntries = 2;
84#endif
85 }
86 else
87 {
88 Indices[0] = MyGlobalElements[i]-1;
89 Indices[1] = MyGlobalElements[i];
90 Indices[2] = MyGlobalElements[i]+1;
91 Values[0] = -1.0;
92 Values[1] = 2.0;
93 Values[2] = -1.0;
94#ifndef NDEBUG
95 NumEntries = 3;
96#endif
97 }
98
99 assert(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
100 // Put in the diagonal entry
101 // assert(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])==0);
102 }
103
104 // Finish up
105 assert(A.FillComplete()==0);
106
107
108 delete[] MyGlobalElements;
109 delete[] Values;
110 delete[] Indices;
111 return 0;
112}
113
114int CreateTridiPlus(Epetra_CrsMatrix& A)
115{
116
117 Epetra_Map Map = A.RowMap();
118 int NumMyElements = Map.NumMyElements();
119 int NumGlobalElements = Map.NumGlobalElements();
120
121 int * MyGlobalElements = new int[NumMyElements];
122 Map.MyGlobalElements(MyGlobalElements);
123
124 // Add rows one-at-a-time
125 // Need some vectors to help
126 // Off diagonal Values will always be -1
127
128
129 double *Values = new double[3];
130 int *Indices = new int[3];
131#ifndef NDEBUG
132 int NumEntries;
133#endif
134 for (int i=0; i<NumMyElements; i++)
135 {
136 if (MyGlobalElements[i]==0)
137 {
138 Indices[0] = 0;
139 Indices[1] = 1;
140 Indices[2] = NumGlobalElements-1;
141 Values[0] = 2.0;
142 Values[1] = -1.0;
143 Values[2] = -0.5;
144#ifndef NDEBUG
145 NumEntries = 3;
146#endif
147 }
148 else if (MyGlobalElements[i] == NumGlobalElements-1)
149 {
150 Indices[0] = NumGlobalElements-1;
151 Indices[1] = NumGlobalElements-2;
152 Indices[2] = 0;
153 Values[0] = 2.0;
154 Values[1] = -1.0;
155 Values[2] = -0.5;
156#ifndef NDEBUG
157 NumEntries = 3;
158#endif
159 }
160 else
161 {
162 Indices[0] = MyGlobalElements[i]-1;
163 Indices[1] = MyGlobalElements[i];
164 Indices[2] = MyGlobalElements[i]+1;
165 Values[0] = -1.0;
166 Values[1] = 2.0;
167 Values[2] = -1.0;
168#ifndef NDEBUG
169 NumEntries = 3;
170#endif
171 }
172
173 assert(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
174 // Put in the diagonal entry
175 // assert(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])==0);
176 }
177
178 // Finish up
179 assert(A.FillComplete()==0);
180
181
182 delete[] MyGlobalElements;
183 delete[] Values;
184 delete[] Indices;
185 return 0;
186}
187
int CreateTridiPlus(Epetra_CrsMatrix &A)
int CreateTridi(Epetra_CrsMatrix &A)
Definition: CreateTridi.cpp:44