Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Euclid_apply.c
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) 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// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Euclid_dh.h"
44#include "Mat_dh.h"
45#include "Factor_dh.h"
46#include "Parser_dh.h"
47#include "TimeLog_dh.h"
48#include "SubdomainGraph_dh.h"
49
50
51static void scale_rhs_private (Euclid_dh ctx, double *rhs);
52static void permute_vec_n2o_private (Euclid_dh ctx, double *xIN,
53 double *xOUT);
54static void permute_vec_o2n_private (Euclid_dh ctx, double *xIN,
55 double *xOUT);
56
57#undef __FUNC__
58#define __FUNC__ "Euclid_dhApply"
59void
60Euclid_dhApply (Euclid_dh ctx, double *rhs, double *lhs)
61{
62 START_FUNC_DH double *rhs_, *lhs_;
63 double t1, t2;
64
65 t1 = MPI_Wtime ();
66
67 /* default settings; for everything except PILU */
68 ctx->from = 0;
69 ctx->to = ctx->m;
70
71 /* case 1: no preconditioning */
72 if (!strcmp (ctx->algo_ilu, "none") || !strcmp (ctx->algo_par, "none"))
73 {
74 int i, m = ctx->m;
75 for (i = 0; i < m; ++i)
76 lhs[i] = rhs[i];
77 goto END_OF_FUNCTION;
78 }
79
80 /*----------------------------------------------------------------
81 * permute and scale rhs vector
82 *----------------------------------------------------------------*/
83 /* permute rhs vector */
84 if (ctx->sg != NULL)
85 {
86
87 permute_vec_n2o_private (ctx, rhs, lhs);
89 rhs_ = lhs;
90 lhs_ = ctx->work2;
91 }
92 else
93 {
94 rhs_ = rhs;
95 lhs_ = lhs;
96 }
97
98 /* scale rhs vector */
99 if (ctx->isScaled)
100 {
101
102 scale_rhs_private (ctx, rhs_);
104 }
105
106 /* note: rhs_ is permuted, scaled; the input, "rhs" vector has
107 not been disturbed.
108 */
109
110 /*----------------------------------------------------------------
111 * big switch to choose the appropriate triangular solve
112 *----------------------------------------------------------------*/
113
114 /* sequential and mpi block jacobi cases */
115 if (np_dh == 1 || !strcmp (ctx->algo_par, "bj"))
116 {
117 Factor_dhSolveSeq (rhs_, lhs_, ctx);
119 }
120
121
122 /* pilu case */
123 else
124 {
125 Factor_dhSolve (rhs_, lhs_, ctx);
127 }
128
129 /*----------------------------------------------------------------
130 * unpermute lhs vector
131 * (note: don't need to unscale, because we were clever)
132 *----------------------------------------------------------------*/
133 if (ctx->sg != NULL)
134 {
135 permute_vec_o2n_private (ctx, lhs_, lhs);
137 }
138
139END_OF_FUNCTION:;
140
141 t2 = MPI_Wtime ();
142 /* collective timing for triangular solves */
143 ctx->timing[TRI_SOLVE_T] += (t2 - t1);
144
145 /* collective timing for setup+krylov+triSolves
146 (intent is to time linear solve, but this is
147 at best probelematical!)
148 */
150
151 /* total triangular solve count */
152 ctx->its += 1;
153 ctx->itsTotal += 1;
154
156
157
158#undef __FUNC__
159#define __FUNC__ "scale_rhs_private"
160void
161scale_rhs_private (Euclid_dh ctx, double *rhs)
162{
163 START_FUNC_DH int i, m = ctx->m;
164 REAL_DH *scale = ctx->scale;
165
166 /* if matrix was scaled, must scale the rhs */
167 if (scale != NULL)
168 {
169 for (i = 0; i < m; ++i)
170 {
171 rhs[i] *= scale[i];
172 }
173 }
175
176
177#undef __FUNC__
178#define __FUNC__ "permute_vec_o2n_private"
179void
180permute_vec_o2n_private (Euclid_dh ctx, double *xIN, double *xOUT)
181{
182 START_FUNC_DH int i, m = ctx->m;
183 int *o2n = ctx->sg->o2n_col;
184 for (i = 0; i < m; ++i)
185 xOUT[i] = xIN[o2n[i]];
187
188
189#undef __FUNC__
190#define __FUNC__ "permute_vec_n2o_private"
191void
192permute_vec_n2o_private (Euclid_dh ctx, double *xIN, double *xOUT)
193{
194 START_FUNC_DH int i, m = ctx->m;
195 int *n2o = ctx->sg->n2o_row;
196 for (i = 0; i < m; ++i)
197 xOUT[i] = xIN[n2o[i]];
void Euclid_dhApply(Euclid_dh ctx, double *rhs, double *lhs)
Definition: Euclid_apply.c:60
static void permute_vec_o2n_private(Euclid_dh ctx, double *xIN, double *xOUT)
Definition: Euclid_apply.c:180
static void permute_vec_n2o_private(Euclid_dh ctx, double *xIN, double *xOUT)
Definition: Euclid_apply.c:192
static void scale_rhs_private(Euclid_dh ctx, double *rhs)
Definition: Euclid_apply.c:161
@ TRI_SOLVE_T
Definition: Euclid_dh.h:111
@ SOLVE_START_T
Definition: Euclid_dh.h:110
@ TOTAL_SOLVE_TEMP_T
Definition: Euclid_dh.h:118
void Factor_dhSolveSeq(double *rhs, double *lhs, Euclid_dh ctx)
Definition: Factor_dh.c:1257
void Factor_dhSolve(double *rhs, double *lhs, Euclid_dh ctx)
Definition: Factor_dh.c:796
int np_dh
Definition: globalObjects.c:62
#define REAL_DH
Definition: euclid_common.h:53
#define START_FUNC_DH
Definition: macros_dh.h:181
#define CHECK_V_ERROR
Definition: macros_dh.h:138
#define END_FUNC_DH
Definition: macros_dh.h:187
SubdomainGraph_dh sg
Definition: Euclid_dh.h:153
double * work2
Definition: Euclid_dh.h:160
char algo_par[MAX_OPT_LEN]
Definition: Euclid_dh.h:164
double timing[TIMING_BINS]
Definition: Euclid_dh.h:189
REAL_DH * scale
Definition: Euclid_dh.h:155
char algo_ilu[MAX_OPT_LEN]
Definition: Euclid_dh.h:165