Sacado Package Browser (Single Doxygen Collection): Sacado_PCE

00001 // $Id: Sacado_PCE_WorkspaceImp.hpp,v 1.1 2007/11/14 00:18:19 etphipp Exp $ 
00002 // $Source: /space/CVS/Trilinos/packages/sacado/src/pce/Sacado_PCE_WorkspaceImp.hpp,v $ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 template <typename BasisT>
00033 Sacado::PCE::Workspace<BasisT>::
00034 Workspace(unsigned int sz_) :
00035   sz(sz_),
00036   A(2*sz,2*sz),
00037   b(2*sz,2),
00038   piv(2*sz),
00039   Cijk(sz-1),
00040   lapack()
00041 {
00042 }
00043 
00044 template <typename BasisT>
00045 void
00046 Sacado::PCE::Workspace<BasisT>::
00047 resize(unsigned int sz_)
00048 {
00049   if (sz_ != sz) {
00050     sz = sz_;
00051     A.shape(2*sz,2*sz);
00052     b.shape(2*sz,2);
00053     piv.resize(2*sz);
00054     Cijk.resize(sz-1);
00055   }
00056 }
00057 
00058 extern "C" {
00059   double dlange_(char*, int*, int*, double*, int*, double*);
00060 }
00061 
00062 template <typename BasisT>
00063 typename Sacado::PCE::Workspace<BasisT>::ordinal_type
00064 Sacado::PCE::Workspace<BasisT>::
00065 solve(typename Sacado::PCE::Workspace<BasisT>::ordinal_type s,
00066       typename Sacado::PCE::Workspace<BasisT>::ordinal_type nrhs)
00067 {
00068   ordinal_type info;
00069 //   lapack.GESV(s, nrhs, A.values(), A.numRows(), &(piv[0]), b.values(), 
00070 //        b.numRows(), &info);
00071   lapack.GETRF(s, s, A.values(), A.numRows(), &(piv[0]), &info);
00072   value_type norm, rcond;
00073   std::vector<ordinal_type> iwork(4*s);
00074   std::vector<value_type> work(4*s);
00075   norm = 1.0;
00076   ordinal_type n = A.numRows();
00077   char t = '1';
00078   norm = dlange_(&t, &s, &s, A.values(), &n, &work[0]);
00079   lapack.GECON('1', s, A.values(), A.numRows(), norm, &rcond, &work[0], 
00080          &iwork[0], &info);
00081   std::cout << "condition number = " << 1.0/rcond << std::endl;
00082   lapack.GETRS('N', s, nrhs, A.values(), A.numRows(), &(piv[0]), b.values(), 
00083          b.numRows(), &info);
00084   return info;
00085 }