Question about how to do something with Pooma 2

[Dave Nystrom]

I have a question which I hope someone can help me with.  I posed this
question to Scott Haney 2-3 years ago in the context of Pooma 1 and he said
that what I wanted to do would be trivial to do with Pooma 2.  Below, I have
included a couple of files, PCG_Ctrl.hh and PCG_Ctrl.cc, which I wrote awhile 
back to interface to a Fortran 77 linear solver package.  One of the key
methods of the PCG_Ctrl class is the pcg_fe method which has the following
signature.

    void pcg_fe( Mat1<T>& x, Mat1<T>& b,
                 SP< PCG_MatVec<T> >& pcg_matvec,
                 SP< PCG_PreCond<T> >& pcg_precond );

The Mat1 class is a data container class that Geoff Furnish wrote and it is
templated on numeric precision type - such as double.  It is a one
dimensional array class and it had a constructor which he called an aliasing
constructor.  The aliasing constructor took a reference to a point in memory
and an integer value that told how many elements of data you were dealing
with.  For a distributed field quantity, the Mat1s were pointing to the first 
element of data that was local to the processor it was on.  The pcg package
has a large work array called fwk which is also a Mat1 but has space
allocated for lots of chunks of data equal to the amount of data which a
field has on a given processor as well as other needed workspace.  For the
way in which I was using the pcg package, it required the user to supply his
own function to perform a matrix-vector multiply operation and a
preconditioning operation.  To do this, the pcg package returns to the
calling function and depending on the value of a flag, it gives you the
offset into the 1D fwk array from a fortran perspective that identifies the
location of the first element of the vector that it wants you to perform some 
operation on and the offset into the 1D fwk array that identifies the first
element where the result vector goes.

I would like to change the signature of the pcg_fe method and other methods
to use Pooma 2 fields i.e. I would like the signature of pcg_fe to become:

    void pcg_fe( Field<...>& x, Field<...>& b,
                 SP< PCG_MatVec<T> >& pcg_matvec,
                 SP< PCG_PreCond<T> >& pcg_precond );

and, I would like fwk to be a 1D Pooma Field and be able to set up views into 
the fwk array that are similar to the way I am doing with the Mat1 class.
Scott thought it would be easy to do this with the Pooma 2 View capability.
I've discussed this with John Hall and he is not sure how to do this off the
top of his head.  One thing we need to do is make sure the guard cell data is 
not in the Field I am using with pcg because pcg does not know anything about 
guard cell data.  However, I will need guard cell data in the matvec function 
and also possibly in the preconditioner functions.  This probably means that
I need at some point to copy from a Field with guard cells to one without and 
vice versa.

Can anyone offer suggestions about how to do what I want using Pooma 2?  Is
my explanation of what I want to do clear enough?  Looking at the pcg_fe
function might also help to clarify things.

Thanks for any advice you can offer.  This will be my first real attempt to
use Pooma 2.

-- 
Dave Nystrom                    email: address@hidden
LANL X-3                        phone: 505-667-7913     fax: 505-665-3046

//----------------------------------*-C++-*----------------------------------//
// PCG_Ctrl.hh
// Dave Nystrom
// Mon Jan 13 17:40:28 1997
//---------------------------------------------------------------------------//
// @> 
//---------------------------------------------------------------------------//

#ifndef __fldeqns_PCG_Ctrl_hh__
#define __fldeqns_PCG_Ctrl_hh__

#include "Mat.hh"
#include "SP.hh"

#include "fldeqns/pcg_DB.hh"
#include "fldeqns/PCG_MatVec.hh"
#include "fldeqns/PCG_PreCond.hh"

//===========================================================================//
// class PCG_Ctrl - 

// 
//===========================================================================//

template<class T>
class PCG_Ctrl : public pcg_DB {
    Mat1<int>            iparm;
    Mat1<T>              fparm;
    Mat1<int>            iwk;
    Mat1<T>              fwk;
    Mat1<T>              xex;
    int                  nru;
    int                  ijob, ireq, iva, ivql, ivqr, ier;
    int                  itmeth, imatvec;

  public:
// Constructor.
    PCG_Ctrl( pcg_DB& _pcg_db, int _nru );

// Destructor.
    ~PCG_Ctrl();

// Public Methods
  public:
    void pcg_fe( Mat1<T>& x, Mat1<T>& b,
                 SP< PCG_MatVec<T> >& pcg_matvec,
                 SP< PCG_PreCond<T> >& pcg_precond );

// Private Methods
  private:
    void set_default();
    void set_nwi();
    void set_nwf();
    void it_method( Mat1<T>& x, Mat1<T>& b, Mat1<T>& xex );
    void print_params();
};

#endif                          // __fldeqns_PCG_Ctrl_hh__

//---------------------------------------------------------------------------//
//                              end of fldeqns/PCG_Ctrl.hh
//---------------------------------------------------------------------------//
//----------------------------------*-C++-*----------------------------------//
// PCG_Ctrl.cc
// Dave Nystrom
// Mon Jan 13 17:40:29 1997
//---------------------------------------------------------------------------//
// @> 
//---------------------------------------------------------------------------//

#include "fldeqns/PCG_Ctrl.hh"
#include "fldeqns/PCG_Subroutines.hh"

//---------------------------------------------------------------------------//
// Destructor.
//---------------------------------------------------------------------------//

template<class T>
PCG_Ctrl<T>::~PCG_Ctrl() {}

//---------------------------------------------------------------------------//
// Constructor.
//---------------------------------------------------------------------------//

template<class T>
PCG_Ctrl<T>::PCG_Ctrl( pcg_DB& pcg_db, int _nru )
    : pcg_DB(pcg_db),
      iparm(Bounds(1,50)), fparm(Bounds(1,30)),
      iwk(pcg_db.nwi), fwk(pcg_db.nwf),
      xex(1), nru(_nru)
{
// Test print.
    cout << "Entering PCG_Ctrl::PCG_Ctrl." << endl;

// Initialize some stuff from pcg_DB.
    itmeth = pcg_db.itmeth;

// Compute required size of iwk and fwk and resize them.
    set_nwi();
    set_nwf();

    iwk.redim( nwi );
    fwk.redim( nwf );

// Initialize iparm and fparm arrays.
    set_default();

    iparm(pcg::NOUT)   = nout;
    iparm(pcg::LEVOUT) = levout;
    iparm(pcg::NRU)    = nru;
    iparm(pcg::ITSMAX) = itsmax;
    iparm(pcg::MALLOC) = malloc;
    iparm(pcg::NWI)    = nwi;
    iparm(pcg::NWF)    = nwf;
    iparm(pcg::NTEST)  = ntest;
    iparm(pcg::IQSIDE) = iqside;
    iparm(pcg::IUINIT) = iuinit;
    iparm(pcg::NEEDRC) = needrc;
    iparm(pcg::NS1)    = ns1;
    iparm(pcg::NS2)    = ns2;
    iparm(pcg::ICKSTG) = ickstg;
    iparm(pcg::IUEXAC) = iuexac;
    iparm(pcg::IDOT)   = idot;
    iparm(pcg::ISTATS) = istats;

    fparm(pcg::CTIMER) = ctimer;
    fparm(pcg::RTIMER) = rtimer;
    fparm(pcg::FLOPSR) = flopsr;
    fparm(pcg::ZETA)   = zeta;
    fparm(pcg::ALPHA)  = alpha;

// Print control parameters.
    print_params();

    cout << "Done with PCG_Ctrl::PCG_Ctrl." << endl << flush;
}

//---------------------------------------------------------------------------//
// Main controller method.
//---------------------------------------------------------------------------//

template<class T>
void PCG_Ctrl<T>::pcg_fe( Mat1<T>& x, Mat1<T>& b,
                          SP< PCG_MatVec<T> >& pcg_matvec,
                          SP< PCG_PreCond<T> >& pcg_precond )
{
// Test print.
    cout << "Entering PCG_Ctrl::pcg_fe." << endl;

// Initialize ijob.
    ijob = pcg::JINIT;

// Call an iterative method.
    int done=0;
    while (!done) {
        it_method( x, b, xex );

        ijob = pcg::JRUN;

        if( ireq == pcg::JTERM ) {
            done = 1;
        }
        else if( ireq == pcg::JAV ) {
            cout << "Preparing for MatVec." << endl << flush;
            Mat1<T> xmatvec(&fwk(ivqr-1),nru);
            Mat1<T> bmatvec(&fwk(iva-1), nru);
            pcg_matvec->MatVec(bmatvec,xmatvec);
            cout << "Done with     MatVec." << endl << flush;
        }
        else if( ireq == pcg::JQLV ) {
            cout << "Preparing for Left_PreCond." << endl << flush;
            Mat1<T> xprecond(&fwk(ivql-1),nru);
            Mat1<T> bprecond(&fwk(iva-1), nru);
            pcg_precond->Left_PreCond(xprecond,bprecond);
            cout << "Done with     Left_PreCond." << endl << flush;
        }
        else if( ireq == pcg::JQRV ) {
            Mat1<T> xprecond(&fwk(ivqr-1),nru);
            Mat1<T> bprecond(&fwk(ivql-1),nru);
            pcg_precond->Right_PreCond(xprecond,bprecond);
        }
        else if( ireq == pcg::JTEST ) {
        }
        else if( ireq == pcg::JATV ) {
        }
        else if( ireq == pcg::JQLTV ) {
        }
        else if( ireq == pcg::JQRTV ) {
        }
    }
}

//---------------------------------------------------------------------------//
// Set default values for PCG iparm and fparm arrays.
//---------------------------------------------------------------------------//

template<class T>
void PCG_Ctrl<T>::set_default()
{
    pcg::xdfalt( &iparm(1), &fparm(1) );
}

//---------------------------------------------------------------------------//
// Call a pcg iterative method.
//---------------------------------------------------------------------------//

template<class T>
void PCG_Ctrl<T>::it_method( Mat1<T>& x, Mat1<T>& b, Mat1<T>& xex )
{
    if( itmeth == pcg::BASIC ) {
        pcg::xbasr( ijob, ireq, &x(0), &xex(0), &b(0), iva, ivql, ivqr,
                    &iwk(0), &fwk(0), &iparm(1), &fparm(1), ier );
        imatvec = pcg::TRUE;
    }
    else if( itmeth == pcg::GMRES ) {
        pcg::xgmrsr( ijob, ireq, &x(0), &xex(0), &b(0), iva, ivql, ivqr,
                     &iwk(0), &fwk(0), &iparm(1), &fparm(1), ier );
        imatvec = pcg::TRUE;
    }
    else {
        throw("Need to choose a valid pcg iterative method.");
    }
}

//---------------------------------------------------------------------------//
// Set nwi.
//---------------------------------------------------------------------------//

template<class T>
void PCG_Ctrl<T>::set_nwi()
{
    if( itmeth == pcg::BASIC    ) nwi = 100;
    if( itmeth == pcg::BCGSTAB  ) nwi = 100;
    if( itmeth == pcg::BCGSTAB2 ) nwi = 100;
    if( itmeth == pcg::BCGSTABL ) nwi = 100;
    if( itmeth == pcg::CGS      ) nwi = 100;
    if( itmeth == pcg::TFQMR    ) nwi = 100;
    if( itmeth == pcg::GMRES    ) nwi = 100;
    if( itmeth == pcg::GMRES_H  ) nwi = 100;
    if( itmeth == pcg::OMIN     ) nwi = 100;
    if( itmeth == pcg::ORES     ) nwi = 100;
    if( itmeth == pcg::IOM      ) nwi = 100;
    if( itmeth == pcg::CG       ) nwi = 100;
    if( itmeth == pcg::BCG      ) nwi = 100;
    if( itmeth == 91            ) nwi = 100;
    if( itmeth == 92            ) nwi = 100;
    if( itmeth == 96            ) nwi = 100;

    if( malloc == pcg::TRUE     ) nwi = 1;

    cout << "------------------------------------------" << endl;
    cout << "Test print from PCG_Ctrl<T>::set_nwi."      << endl;
    cout << "------------------------------------------" << endl;
    cout << "     itmeth = " << itmeth << endl;
    cout << "     nwi    = " << nwi    << endl;
}

//---------------------------------------------------------------------------//
// Set nwf.
//---------------------------------------------------------------------------//

template<class T>
void PCG_Ctrl<T>::set_nwf()
{
    int nrup2   = nru + 2;
    int nwfgenl = 31 + nrup2*2;
    int nwfit   = 0;
    int nwfstat = 0;
    int nwftst  = 0;

    if( itmeth == pcg::BASIC    ) nwfit =  7 + nrup2*5;
    if( itmeth == pcg::BCGSTAB  ) nwfit = 15 + nrup2*13;
    if( itmeth == pcg::BCGSTAB2 ) nwfit = 15 + nrup2*25;
    if( itmeth == pcg::BCGSTABL ) nwfit = 31 + nrup2*(2*ns2+8) + ns2*9
                                        + ns2*ns2;
    if( itmeth == pcg::CGS      ) nwfit = 12 + nrup2*11;
    if( itmeth == pcg::TFQMR    ) nwfit = 16 + nrup2*18;
    if( itmeth == pcg::CG       ) nwfit = 12 + nrup2*5;
    if( itmeth == pcg::BCG      ) nwfit = 12 + nrup2*9;
    if( itmeth == 91            ) nwfit =  7 + nrup2*5;
    if( itmeth == 92            ) nwfit = 15 + nrup2*13;
    if( iqside >= pcg::QRIGHT ) {
        if( itmeth == pcg::GMRES   ) nwfit = 31 + nrup2*(2*ns2+8) + ns2*9
                                           + ns2*ns2;
        if( itmeth == pcg::GMRES_H ) nwfit = 26 + nrup2*(2*ns2+7) + ns2*7
                                           + ns2*ns2;
        if( itmeth == pcg::OMIN    ) nwfit = 13 + nrup2*(3*ns1+5) + ns1;
        if( itmeth == pcg::ORES    ) nwfit = 12 + nrup2*(4*ns1+6) + ns1;
        if( itmeth == pcg::IOM     ) nwfit = 31 + nrup2*(3*ns1+9) + ns1*5;
        if( itmeth == 96           ) nwfit = 31 + nrup2*(2*ns2+8) + ns2*9
                                           + ns2*ns2;
    }        
    else {
        if( itmeth == pcg::GMRES   ) nwfit = 31 + nrup2*(1*ns2+4) + ns2*9
                                           + ns2*ns2;
        if( itmeth == pcg::GMRES_H ) nwfit = 26 + nrup2*(2*ns2+7) + ns2*7
                                           + ns2*ns2;
        if( itmeth == pcg::OMIN    ) nwfit = 13 + nrup2*(2*ns1+3) + ns1;
        if( itmeth == pcg::ORES    ) nwfit = 12 + nrup2*(2*ns1+4) + ns1;
        if( itmeth == pcg::IOM     ) nwfit = 31 + nrup2*(2*ns1+7) + ns1*5;
        if( itmeth == 96           ) nwfit = 31 + nrup2*(1*ns2+4) + ns2*9
                                           + ns2*ns2;
    }

    if( istats == 1 ) nwfstat = 20 + nrup2*4;

    if( ntest != pcg::TST0 && ntest != pcg::TSTDFA ) {
        nwftst = nrup2*2;
    }

    if( malloc == pcg::TRUE ) {
        nwf = 1;
    }
    else {
        nwf = nwfgenl + nwfstat + nwftst + nwfit;
    }

    cout << "------------------------------------------" << endl;
    cout << "Test print from PCG_Ctrl<T>::set_nwf."      << endl;
    cout << "------------------------------------------" << endl;
    cout << "     itmeth = " << itmeth << endl;
    cout << "     nwf    = " << nwf    << endl;

//    nwf = 1000000;
}

//---------------------------------------------------------------------------//
// Print values for PCG iparm and fparm arrays.
//---------------------------------------------------------------------------//

template<class T>
void PCG_Ctrl<T>::print_params()
{
// Test print.
    cout << "Entering PCG_Ctrl::print_params." << endl;

// Revcom level parameters.
    cout << "----------------------------------------------" << endl;
    cout << "Revcom level parameters."                       << endl;
    cout << "----------------------------------------------" << endl;
    cout << "     nout   = " << iparm(pcg::NOUT)   << endl;
    cout << "     levout = " << iparm(pcg::LEVOUT) << endl;
    cout << "     nru    = " << iparm(pcg::NRU)    << endl;
    cout << "     itsmax = " << iparm(pcg::ITSMAX) << endl;
    cout << "     its    = " << iparm(pcg::ITS)    << endl;
    cout << "     malloc = " << iparm(pcg::MALLOC) << endl;
    cout << "     nwi    = " << iparm(pcg::NWI)    << endl;
    cout << "     nwf    = " << iparm(pcg::NWF)    << endl;
    cout << "     nwiusd = " << iparm(pcg::NWIUSD) << endl;
    cout << "     nwfusd = " << iparm(pcg::NWFUSD) << endl;
    cout << "     iptr   = " << iparm(pcg::IPTR)   << endl;
    cout << "     ntest  = " << iparm(pcg::NTEST)  << endl;
    cout << "     iqside = " << iparm(pcg::IQSIDE) << endl;
    cout << "     iuinit = " << iparm(pcg::IUINIT) << endl;
    cout << "     needrc = " << iparm(pcg::NEEDRC) << endl;
    cout << "     ns1    = " << iparm(pcg::NS1)    << endl;
    cout << "     ns2    = " << iparm(pcg::NS2)    << endl;
    cout << "     ickstg = " << iparm(pcg::ICKSTG) << endl;
    cout << "     iuexac = " << iparm(pcg::IUEXAC) << endl;
    cout << "     idot   = " << iparm(pcg::IDOT)   << endl;
    cout << "     istats = " << iparm(pcg::ISTATS) << endl;
    cout << "     itimer = " << iparm(pcg::ITIMER) << endl;
    cout << "     icomm  = " << iparm(pcg::ICOMM)  << endl;
    cout << "     msgmin = " << iparm(pcg::MSGMIN) << endl;
    cout << "     msgmax = " << iparm(pcg::MSGMAX) << endl;
    cout << "     msgtyp = " << iparm(pcg::MSGTYP) << endl;
    cout << "     iclev  = " << iparm(pcg::ICLEV)  << endl;
    cout << " "                                    << endl;
    cout << "     ctimer = " << fparm(pcg::CTIMER) << endl;
    cout << "     rtimer = " << fparm(pcg::RTIMER) << endl;
    cout << "     flopsr = " << fparm(pcg::FLOPSR) << endl;
    cout << "     zeta   = " << fparm(pcg::ZETA)   << endl;
    cout << "     stptst = " << fparm(pcg::STPTST) << endl;
    cout << "     alpha  = " << fparm(pcg::ALPHA)  << endl;
    cout << "     relrsd = " << fparm(pcg::RELRSD) << endl;
    cout << "     relerr = " << fparm(pcg::RELERR) << endl;
}

//---------------------------------------------------------------------------//
//                              end of PCG_Ctrl.cc
//---------------------------------------------------------------------------//