IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Euclid_apply.c
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 
51 static void scale_rhs_private (Euclid_dh ctx, double *rhs);
52 static void permute_vec_n2o_private (Euclid_dh ctx, double *xIN,
53  double *xOUT);
54 static void permute_vec_o2n_private (Euclid_dh ctx, double *xIN,
55  double *xOUT);
56 
57 #undef __FUNC__
58 #define __FUNC__ "Euclid_dhApply"
59 void
60 Euclid_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);
88  CHECK_V_ERROR;
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_);
103  CHECK_V_ERROR;
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);
118  CHECK_V_ERROR;
119  }
120 
121 
122  /* pilu case */
123  else
124  {
125  Factor_dhSolve (rhs_, lhs_, ctx);
126  CHECK_V_ERROR;
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);
136  CHECK_V_ERROR;
137  }
138 
139 END_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  */
149  ctx->timing[TOTAL_SOLVE_TEMP_T] = t2 - ctx->timing[SOLVE_START_T];
150 
151  /* total triangular solve count */
152  ctx->its += 1;
153  ctx->itsTotal += 1;
154 
155 END_FUNC_DH}
156 
157 
158 #undef __FUNC__
159 #define __FUNC__ "scale_rhs_private"
160 void
161 scale_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  }
174 END_FUNC_DH}
175 
176 
177 #undef __FUNC__
178 #define __FUNC__ "permute_vec_o2n_private"
179 void
180 permute_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]];
186 END_FUNC_DH}
187 
188 
189 #undef __FUNC__
190 #define __FUNC__ "permute_vec_n2o_private"
191 void
192 permute_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]];
198 END_FUNC_DH}