Euclid_apply.c

00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Euclid_dh.h"
00031 #include "Mat_dh.h"
00032 #include "Factor_dh.h"
00033 #include "Parser_dh.h"
00034 #include "TimeLog_dh.h"
00035 #include "SubdomainGraph_dh.h"
00036 
00037 
00038 static void scale_rhs_private (Euclid_dh ctx, double *rhs);
00039 static void permute_vec_n2o_private (Euclid_dh ctx, double *xIN,
00040                      double *xOUT);
00041 static void permute_vec_o2n_private (Euclid_dh ctx, double *xIN,
00042                      double *xOUT);
00043 
00044 #undef __FUNC__
00045 #define __FUNC__ "Euclid_dhApply"
00046 void
00047 Euclid_dhApply (Euclid_dh ctx, double *rhs, double *lhs)
00048 {
00049   START_FUNC_DH double *rhs_, *lhs_;
00050   double t1, t2;
00051 
00052   t1 = MPI_Wtime ();
00053 
00054   /* default settings; for everything except PILU */
00055   ctx->from = 0;
00056   ctx->to = ctx->m;
00057 
00058   /* case 1: no preconditioning */
00059   if (!strcmp (ctx->algo_ilu, "none") || !strcmp (ctx->algo_par, "none"))
00060     {
00061       int i, m = ctx->m;
00062       for (i = 0; i < m; ++i)
00063     lhs[i] = rhs[i];
00064       goto END_OF_FUNCTION;
00065     }
00066 
00067   /*----------------------------------------------------------------
00068    * permute and scale rhs vector
00069    *----------------------------------------------------------------*/
00070   /* permute rhs vector */
00071   if (ctx->sg != NULL)
00072     {
00073 
00074       permute_vec_n2o_private (ctx, rhs, lhs);
00075       CHECK_V_ERROR;
00076       rhs_ = lhs;
00077       lhs_ = ctx->work2;
00078     }
00079   else
00080     {
00081       rhs_ = rhs;
00082       lhs_ = lhs;
00083     }
00084 
00085   /* scale rhs vector */
00086   if (ctx->isScaled)
00087     {
00088 
00089       scale_rhs_private (ctx, rhs_);
00090       CHECK_V_ERROR;
00091     }
00092 
00093   /* note: rhs_ is permuted, scaled; the input, "rhs" vector has
00094      not been disturbed.
00095    */
00096 
00097   /*----------------------------------------------------------------
00098    * big switch to choose the appropriate triangular solve
00099    *----------------------------------------------------------------*/
00100 
00101   /* sequential and mpi block jacobi cases */
00102   if (np_dh == 1 || !strcmp (ctx->algo_par, "bj"))
00103     {
00104       Factor_dhSolveSeq (rhs_, lhs_, ctx);
00105       CHECK_V_ERROR;
00106     }
00107 
00108 
00109   /* pilu case */
00110   else
00111     {
00112       Factor_dhSolve (rhs_, lhs_, ctx);
00113       CHECK_V_ERROR;
00114     }
00115 
00116   /*----------------------------------------------------------------
00117    * unpermute lhs vector
00118    * (note: don't need to unscale, because we were clever)
00119    *----------------------------------------------------------------*/
00120   if (ctx->sg != NULL)
00121     {
00122       permute_vec_o2n_private (ctx, lhs_, lhs);
00123       CHECK_V_ERROR;
00124     }
00125 
00126 END_OF_FUNCTION:;
00127 
00128   t2 = MPI_Wtime ();
00129   /* collective timing for triangular solves */
00130   ctx->timing[TRI_SOLVE_T] += (t2 - t1);
00131 
00132   /* collective timing for setup+krylov+triSolves
00133      (intent is to time linear solve, but this is
00134      at best probelematical!)
00135    */
00136   ctx->timing[TOTAL_SOLVE_TEMP_T] = t2 - ctx->timing[SOLVE_START_T];
00137 
00138   /* total triangular solve count */
00139   ctx->its += 1;
00140   ctx->itsTotal += 1;
00141 
00142 END_FUNC_DH}
00143 
00144 
00145 #undef __FUNC__
00146 #define __FUNC__ "scale_rhs_private"
00147 void
00148 scale_rhs_private (Euclid_dh ctx, double *rhs)
00149 {
00150   START_FUNC_DH int i, m = ctx->m;
00151   REAL_DH *scale = ctx->scale;
00152 
00153   /* if matrix was scaled, must scale the rhs */
00154   if (scale != NULL)
00155     {
00156       for (i = 0; i < m; ++i)
00157     {
00158       rhs[i] *= scale[i];
00159     }
00160     }
00161 END_FUNC_DH}
00162 
00163 
00164 #undef __FUNC__
00165 #define __FUNC__ "permute_vec_o2n_private"
00166 void
00167 permute_vec_o2n_private (Euclid_dh ctx, double *xIN, double *xOUT)
00168 {
00169   START_FUNC_DH int i, m = ctx->m;
00170   int *o2n = ctx->sg->o2n_col;
00171   for (i = 0; i < m; ++i)
00172     xOUT[i] = xIN[o2n[i]];
00173 END_FUNC_DH}
00174 
00175 
00176 #undef __FUNC__
00177 #define __FUNC__ "permute_vec_n2o_private"
00178 void
00179 permute_vec_n2o_private (Euclid_dh ctx, double *xIN, double *xOUT)
00180 {
00181   START_FUNC_DH int i, m = ctx->m;
00182   int *n2o = ctx->sg->n2o_row;
00183   for (i = 0; i < m; ++i)
00184     xOUT[i] = xIN[n2o[i]];
00185 END_FUNC_DH}

Generated on Wed May 12 21:30:17 2010 for IFPACK by  doxygen 1.4.7