Commit 527b1490 authored by Francois Gygi's avatar Francois Gygi
Browse files

removed IntegerMatrix class

git-svn-id: http://qboxcode.org/svn/qb/trunk@1098 cba15fb0-1239-40c8-b417-11db7ca47a34
parent cd100b71
......@@ -15,7 +15,6 @@
// Matrix.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Matrix.C,v 1.22 2010-05-10 20:50:37 fgygi Exp $
#include <cassert>
#include <vector>
......@@ -71,9 +70,10 @@ using namespace std;
#define pdlacp3 pdlacp3_
#define pdgetrf pdgetrf_
#define pdgetri pdgetri_
#define pilapiv pilapiv_
#define pdlapiv pdlapiv_
#define pzlapiv pzlapiv_
#define pdlapv2 pdlapv2_
#define pzlapv2 pzlapv2_
#define dscal dscal_
#define zscal zscal_
......@@ -248,6 +248,7 @@ extern "C"
void pdgetri(const int* n, double* val,
const int* ia, const int* ja, int* desca, int* ipiv,
double* work, int* lwork, int* iwork, int* liwork, int* info);
void pdlapiv(const char* direc, const char* rowcol, const char* pivroc,
const int* m, const int* n, double *a, const int* ia,
const int* ja, const int* desca, int* ipiv, const int* ip,
......@@ -256,6 +257,14 @@ extern "C"
const int* m, const int* n, complex<double> *a, const int* ia,
const int* ja, const int* desca, int* ipiv, const int* ip,
const int* jp, const int* descp, int* iwork);
void pdlapv2(const char* direc, const char *rowcol,
const int* m, const int *n, double *val,
const int *ia, const int *ja, const int* desca,
int *ipiv, const int *ip, const int *jp, const int *descp);
void pzlapv2(const char* direc, const char *rowcol,
const int* m, const int *n, complex<double> *val,
const int *ia, const int *ja, const int* desca,
int *ipiv, const int *ip, const int *jp, const int *descp);
#endif
// BLAS1
......@@ -404,76 +413,6 @@ ComplexMatrix::ComplexMatrix(const DoubleMatrix& rhs) : ctxt_(rhs.context()),
val = (complex<double>*) rhs.cvalptr();
}
////////////////////////////////////////////////////////////////////////////////
void IntegerMatrix::init_size(int m, int n, int mb, int nb)
{
assert(m>=0);
assert(n>=0);
assert(mb>=0);
assert(nb>=0);
m_ = m;
n_ = n;
#ifdef SCALAPACK
mb_ = mb;
nb_ = nb;
#else
mb_ = m;
nb_ = n;
#endif
if ( mb_ == 0 ) mb_ = 1;
if ( nb_ == 0 ) nb_ = 1;
ictxt_ = ctxt_.ictxt();
nprow_ = ctxt_.nprow();
npcol_ = ctxt_.npcol();
myrow_ = ctxt_.myrow();
mycol_ = ctxt_.mycol();
active_ = myrow_ >= 0;
int isrcproc=0;
mloc_ = 0;
nloc_ = 0;
if ( m_ != 0 )
mloc_ = numroc(&m_,&mb_,&myrow_,&isrcproc,&nprow_);
if ( n_ != 0 )
nloc_ = numroc(&n_,&nb_,&mycol_,&isrcproc,&npcol_);
size_ = mloc_ * nloc_;
// set leading dimension of val array to mloc_;
lld_ = mloc_;
if ( lld_ == 0 ) lld_ = 1;
// total and local number of blocks
mblocks_ = 0;
nblocks_ = 0;
m_incomplete_ = false;
n_incomplete_ = false;
if ( active_ && mb_ > 0 && nb_ > 0 )
{
mblocks_ = ( mloc_ + mb_ - 1 ) / mb_;
nblocks_ = ( nloc_ + nb_ - 1 ) / nb_;
m_incomplete_ = mloc_ % mb_ != 0;
n_incomplete_ = nloc_ % nb_ != 0;
}
if ( active_ )
{
desc_[0] = 1;
}
else
{
desc_[0] = -1;
}
desc_[1] = ictxt_;
desc_[2] = m_;
desc_[3] = n_;
desc_[4] = mb_;
desc_[5] = nb_;
desc_[6] = 0;
desc_[7] = 0;
desc_[8] = lld_;
}
////////////////////////////////////////////////////////////////////////////////
void DoubleMatrix::init_size(int m, int n, int mb, int nb)
{
......@@ -540,7 +479,6 @@ void DoubleMatrix::init_size(int m, int n, int mb, int nb)
desc_[6] = 0;
desc_[7] = 0;
desc_[8] = lld_;
}
////////////////////////////////////////////////////////////////////////////////
......@@ -610,14 +548,6 @@ void ComplexMatrix::init_size(int m, int n, int mb, int nb)
desc_[6] = 0;
desc_[7] = 0;
desc_[8] = lld_;
}
////////////////////////////////////////////////////////////////////////////////
void IntegerMatrix::clear(void)
{
assert(val!=0||size_==0);
memset(val,0,size_*sizeof(int));
}
////////////////////////////////////////////////////////////////////////////////
......@@ -1018,57 +948,6 @@ void ComplexMatrix::axpy(double alpha, const ComplexMatrix &x)
daxpy(&len, &alpha, (double*) x.val, &ione, (double*) val, &ione);
}
////////////////////////////////////////////////////////////////////////////////
// real getsub: *this = sub(A)
// copy submatrix A(ia:ia+m, ja:ja+n) into *this;
// *this and A may live in different contexts
void IntegerMatrix::getsub(const IntegerMatrix &a,
int m, int n, int ia, int ja)
{
#if SCALAPACK
int iap=ia+1;
int jap=ja+1;
assert(n<=n_);
assert(n<=a.n());
assert(m<=m_);
assert(m<=a.m());
int ione = 1;
int gictxt;
Cblacs_get( 0, 0, &gictxt );
pigemr2d(&m,&n,a.val,&iap,&jap,a.desc_,val,&ione,&ione,desc_,&gictxt);
#else
for ( int j = 0; j < n; j++ )
for ( int i = 0; i < m; i++ )
val[i+j*m_] = a.val[(i+ia) + (j+ja)*a.m()];
#endif
}
////////////////////////////////////////////////////////////////////////////////
// real getsub: *this = sub(A)
// copy submatrix A(ia:ia+m, ja:ja+n) into *this(idest:idest+m,jdest:jdest+n)
// *this and A may live in different contexts
void IntegerMatrix::getsub(const IntegerMatrix &a,
int m, int n, int isrc, int jsrc, int idest, int jdest)
{
#if SCALAPACK
int iap=isrc+1;
int jap=jsrc+1;
int idp=idest+1;
int jdp=jdest+1;
assert(n<=n_);
assert(n<=a.n());
assert(m<=m_);
assert(m<=a.m());
int gictxt;
Cblacs_get( 0, 0, &gictxt );
pigemr2d(&m,&n,a.val,&iap,&jap,a.desc_,val,&idp,&jdp,desc_,&gictxt);
#else
for ( int j = 0; j < n; j++ )
for ( int i = 0; i < m; i++ )
val[(idest+i)+(jdest+j)*m_] = a.val[(i+isrc) + (j+jsrc)*a.m()];
#endif
}
////////////////////////////////////////////////////////////////////////////////
// real getsub: *this = sub(A)
// copy submatrix A(ia:ia+m, ja:ja+n) into *this;
......@@ -3097,90 +2976,134 @@ void ComplexMatrix::heevd(char uplo, valarray<double>& w)
#if SCALAPACK
////////////////////////////////////////////////////////////////////////////////
void DoubleMatrix::lapiv(char direc, char rowcol, IntegerMatrix &permutation)
{
// verify alignment
if ( rowcol=='C' || rowcol=='c' )
{
assert( permutation.n()==n_ );
assert( permutation.nb()==nb_ );
assert( permutation.m()==1 );
assert( permutation.mb()==1 );
void DoubleMatrix::lapiv(char direc, char rowcol, const int *ipiv)
{
// Perform a permutation of rows or columns of *this
//
// Permutation of rows: (rowcol='R' or 'r')
// the array ipiv is distributed over a process column
// and is replicated on all process columns
// ipiv has size mloc and contains the local values of the permutation
//
// Permutation of columns: (rowcol='C' or 'c')
// the array ipiv is distributed over a process row
// and is replicated on all process rows
// ipiv has size nloc and contains the local values of the permutation
const bool rowcol_r = ( rowcol=='R' || rowcol=='r' );
const bool rowcol_c = ( rowcol=='C' || rowcol=='c' );
assert(rowcol_r || rowcol_c);
// ipivtmp: extended permutation array for use in pdlapv2
// (see scalapack documentation)
vector<int> ipivtmp;
// descriptor of the ipivtmp distributed vector
int desc_ip[9];
if ( rowcol_r )
{
// permuting rows
ipivtmp.resize(mloc_ + mb_);
for ( int i = 0; i < mloc_; i++)
ipivtmp[i] = ipiv[i];
// initialize descriptor: ipivtmp is (mx1)
desc_ip[0] = 1; // dtype
desc_ip[1] = ictxt_;// ctxt
desc_ip[2] = m_+mb_*nprow_; // m (see details in pldapv2.f)
desc_ip[3] = 1; // n
desc_ip[4] = mb_; // mb
desc_ip[5] = 1; // nb
desc_ip[6] = 0; // rsrc
desc_ip[7] = 0; // csrc
desc_ip[8] = mloc_; // lld
}
else if ( rowcol=='R' || rowcol=='R' )
else
{
assert( permutation.m()==m_ );
assert( permutation.mb()==mb_ );
assert( permutation.n()==1 );
assert( permutation.nb()==1 );
// permuting columns
ipivtmp.resize(nloc_ + nb_);
for ( int i = 0; i < nloc_; i++)
ipivtmp[i] = ipiv[i];
// initialize descriptor: ipivtmp is (1xn)
desc_ip[0] = 1; // dtype
desc_ip[1] = ictxt_;// ctxt
desc_ip[2] = 1; // m
desc_ip[3] = n_+nb_*npcol_; // n (see details in pdlapv2.f)
desc_ip[4] = 1; // mb
desc_ip[5] = nb_; // nb
desc_ip[6] = 0; // rsrc
desc_ip[7] = 0; // csrc
desc_ip[8] = 1; // lld
}
char pivroc=rowcol;
int one = 1;
vector<int> iwork(mb_+nb_);
// apply permutation
pdlapiv(&direc, &rowcol, &pivroc, &m_, &n_, val, &one, &one, desc_,
&permutation[0], &one, &one, permutation.desc(), &iwork[0] );
}
////////////////////////////////////////////////////////////////////////////////
void ComplexMatrix::lapiv(char direc, char rowcol, IntegerMatrix &permutation)
{
// verify alignment
if ( rowcol=='C' || rowcol=='c' )
{
assert( permutation.n()==n_ );
assert( permutation.nb()==nb_ );
assert( permutation.m()==1 );
assert( permutation.mb()==1 );
int one = 1;
pdlapv2(&direc, &rowcol, &m_, &n_, val, &one, &one, desc_,
&ipivtmp[0], &one, &one, desc_ip);
}
////////////////////////////////////////////////////////////////////////////////
void ComplexMatrix::lapiv(char direc, char rowcol, const int *ipiv)
{
// Perform a permutation of rows or columns of *this
//
// Permutation of rows: (rowcol='R' or 'r')
// the array ipiv is distributed over a process column
// and is replicated on all process columns
// ipiv has size mloc and contains the local values of the permutation
//
// Permutation of columns: (rowcol='C' or 'c')
// the array ipiv is distributed over a process row
// and is replicated on all process rows
// ipiv has size nloc and contains the local values of the permutation
const bool rowcol_r = ( rowcol=='R' || rowcol=='r' );
const bool rowcol_c = ( rowcol=='C' || rowcol=='c' );
assert(rowcol_r || rowcol_c);
// ipivtmp: extended permutation array for use in pdlapv2
// (see scalapack documentation)
vector<int> ipivtmp;
// descriptor of the ipivtmp distributed vector
int desc_ip[9];
if ( rowcol_r )
{
// permuting rows
ipivtmp.resize(mloc_ + mb_);
for ( int i = 0; i < mloc_; i++)
ipivtmp[i] = ipiv[i];
// initialize descriptor: ipivtmp is (mx1)
desc_ip[0] = 1; // dtype
desc_ip[1] = ictxt_;// ctxt
desc_ip[2] = m_+mb_*nprow_; // m (see details in pldapv2.f)
desc_ip[3] = 1; // n
desc_ip[4] = mb_; // mb
desc_ip[5] = 1; // nb
desc_ip[6] = 0; // rsrc
desc_ip[7] = 0; // csrc
desc_ip[8] = mloc_; // lld
}
else if ( rowcol=='R' || rowcol=='R' )
else
{
assert( permutation.m()==m_ );
assert( permutation.mb()==mb_ );
assert( permutation.n()==1 );
assert( permutation.nb()==1 );
// permuting columns
ipivtmp.resize(nloc_ + nb_);
for ( int i = 0; i < nloc_; i++)
ipivtmp[i] = ipiv[i];
// initialize descriptor: ipivtmp is (1x(n+nb))
desc_ip[0] = 1; // dtype
desc_ip[1] = ictxt_;// ctxt
desc_ip[2] = 1; // m
desc_ip[3] = n_+nb_*npcol_; // n (see details in pdlapv2.f)
desc_ip[4] = 1; // mb
desc_ip[5] = nb_; // nb
desc_ip[6] = 0; // rsrc
desc_ip[7] = 0; // csrc
desc_ip[8] = 1; // lld
}
int one = 1;
char pivroc=rowcol;
vector<int> iwork(m_+n_);
// apply permutation
pzlapiv(&direc, &rowcol, &pivroc, &m_, &n_, val, &one, &one, desc_,
permutation.valptr(0), &one, &one, permutation.desc(), &iwork[0] );
pzlapv2(&direc, &rowcol, &m_, &n_, val, &one, &one, desc_,
&ipivtmp[0], &one, &one, desc_ip);
}
#endif
////////////////////////////////////////////////////////////////////////////////
void IntegerMatrix::print(ostream& os) const
{
// Copy blocks of <blocksize> columns and print them on process (0,0)
if ( m_ == 0 || n_ == 0 ) return;
Context ctxtl(1,1);
const int blockmemsize = 32768; // maximum memory size of a block in bytes
// compute maximum block size: must be at least 1
int maxbs = max(1, (int) ((blockmemsize/sizeof(double))/m_));
IntegerMatrix t(ctxtl,m_,maxbs);
int nblocks = n_ / maxbs + ( (n_%maxbs == 0) ? 0 : 1 );
int ia = 0;
int ja = 0;
for ( int jb = 0; jb < nblocks; jb++ )
{
int blocksize = ( (jb+1) * maxbs > n_ ) ? n_ % maxbs : maxbs;
t.getsub(*this,t.m(),blocksize,ia,ja);
ja += blocksize;
if ( t.active() )
{
// this is done only on pe 0
for ( int jj = 0; jj < blocksize; jj++ )
{
for ( int ii = 0; ii < m_; ii++ )
{
os << "(" << ii << "," << jj+jb*maxbs << ")="
<< t.val[ii+t.mloc()*jj] << endl;
}
}
}
}
}
////////////////////////////////////////////////////////////////////////////////
void DoubleMatrix::print(ostream& os) const
{
......@@ -3245,12 +3168,7 @@ void ComplexMatrix::print(ostream& os) const
}
}
}
////////////////////////////////////////////////////////////////////////////////
ostream& operator<<(ostream& os, const IntegerMatrix& a)
{
a.print(os);
return os;
}
////////////////////////////////////////////////////////////////////////////////
ostream& operator<<(ostream& os, const DoubleMatrix& a)
{
......
......@@ -15,7 +15,6 @@
// Matrix.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Matrix.h,v 1.23 2010-05-10 20:49:38 fgygi Exp $
#ifndef MATRIX_H
#define MATRIX_H
......@@ -30,159 +29,6 @@ class Context;
class ComplexMatrix;
class IntegerMatrix
{
private:
Context ctxt_;
int ictxt_;
int lld_; // leading dimension of local matrix
int m_, n_; // size of global matrix
int mb_, nb_; // size of blocks
int mloc_, nloc_; // size of local array
int size_; // size of local array;
int nprow_, npcol_; // number of process rows and cols in context
int myrow_, mycol_; // position of my process in the process grid
int mblocks_, nblocks_; // number of local blocks
int desc_[9];
bool active_;
bool m_incomplete_, n_incomplete_; // this process has an incomplete block
bool reference_; // object was created using the copy constructor
int* val;
public:
int* valptr(int i=0) { return &val[i]; }
const int* cvalptr(int i=0) const { return &val[i]; }
int& operator[] (int i) { return val[i]; }
const int& operator[] (int i) const { return val[i]; }
const Context& context(void) const { return ctxt_; }
int ictxt(void) const { return ictxt_; }
int lld(void) const { return lld_; }
int m(void) const { return m_; } // size of global matrix
int n(void) const { return n_; } // size of global matrix
int mb(void) const { return mb_; } // size of blocks
int nb(void) const { return nb_; } // size of blocks
int mloc(void) const { return mloc_; } // size of local array
int nloc(void) const { return nloc_; } // size of local array
int size(void) const { return size_; } // size of local array
int localsize(void) const { return mloc_*nloc_; } // local size of val
double memsize(void) const { return (double)m_ * (double)n_ * sizeof(double); }
double localmemsize(void) const
{ return (double) mloc_ * (double) nloc_ * sizeof(double); }
const int* desc(void) const { return &desc_[0]; }
// local block size of block (l,m)
int mbs(int l) const
{
// if l is the last block and this process holds an incomplete
// block, then the size is m_%mb_. Otherwise, it is mb_.
return ( (l==mblocks_-1) && ( m_incomplete_ ) ) ? m_%mb_ : mb_;
}
int nbs(int m) const
{
// if m is the last block and this process holds an incomplete
// block, then the size is n_%nb_. Otherwise, it is nb_.
return ( (m==nblocks_-1) && ( n_incomplete_ ) ) ? n_%nb_ : nb_;
}
// number of local blocks
int mblocks(void) const { return mblocks_; }
int nblocks(void) const { return nblocks_; }
// functions to compute local indices from global indices
// index of blocks: element (i,j) is at position (x,y)
// in local block (l,m) of process (pr,pc)
int l(int i) const { return i/(nprow_*mb_); }
int x(int i) const { return i % mb_; }
int pr(int i) const { return (i/mb_) % nprow_; }
int m(int j) const { return j/(npcol_*nb_); }
int y(int j) const { return j % nb_; }
int pc(int j) const { return (j/nb_) % npcol_; }
// global indices:
// (i,j) is the global index of element (x,y) of block (l,m)
int i(int l, int x) const { return (l * nprow_ + myrow_) * mb_ + x; }
int j(int m, int y) const { return (m * npcol_ + mycol_) * nb_ + y; }
int iglobal(int ilocal) const
{ return mb_*(nprow_*(ilocal/mb_)+myrow_)+ilocal%mb_; }
int jglobal(int jlocal) const
{ return nb_*(npcol_*(jlocal/nb_)+mycol_)+jlocal%nb_; }
int iglobal(int iprow, int ilocal) const
{ return mb_*(nprow_*(ilocal/mb_)+iprow)+ilocal%mb_; }
int jglobal(int jpcol, int jlocal) const
{ return nb_*(npcol_*(jlocal/nb_)+jpcol)+jlocal%nb_; }
// store element a(ii,jj) (where ii,jj are global indices)
// in array val:
// int iii = l(ii) * mb_ + x(ii);
// int jjj = m(jj) * nb_ + y(jj);
// val[iii+mloc_*jjj] = a_ij;
// active flag: the matrix has elements on this process
bool active(void) const { return active_; }
void init_size(int m, int n, int mb = MATRIX_DEF_BLOCK_SIZE,
int nb = MATRIX_DEF_BLOCK_SIZE);
void resize(int m, int n, int mb = MATRIX_DEF_BLOCK_SIZE,
int nb = MATRIX_DEF_BLOCK_SIZE)
{
const int old_size = size_;
init_size(m,n,mb,nb);
if ( size_ == old_size ) return;
init_size(m,n,mb,nb);
delete[] val;
val = new int[size_];
clear();
}
void clear(void);
void print(std::ostream& os) const;
explicit IntegerMatrix(const Context& ctxt) : ctxt_(ctxt),
m_(0), n_(0), mb_(0), nb_(0), size_(0), reference_(false), val(0) {}
// Construct an IntegerMatrix of dimensions m,n
explicit IntegerMatrix(const Context& ctxt, int m, int n,
int mb=MATRIX_DEF_BLOCK_SIZE,
int nb=MATRIX_DEF_BLOCK_SIZE) : ctxt_(ctxt),
m_(0), n_(0), mb_(0), nb_(0), size_(0), reference_(false), val(0)
{
resize(m,n,mb,nb);
}
// copy constructor: create a separate copy of rhs
explicit IntegerMatrix(const IntegerMatrix& rhs) : ctxt_(rhs.context()),
reference_(false)
{
init_size(rhs.m(),rhs.n(),rhs.mb(),rhs.nb());
val = new int[size_];
memcpy(val, rhs.val, size_*sizeof(int));
}
// get submatrix A(ia:ia+m,ja:ja+n) of A
void getsub(const IntegerMatrix& a,int m,int n,int ia,int ja);
// get submatrix A(ia:ia+m,ja:ja+n) of A and store in
// this(idest:idest+m,jdest:jdest+n)
void getsub(const IntegerMatrix& a,int m,int n,int isrc,int jsrc,
int idest, int jdest);
// destructor
~IntegerMatrix(void)
{
if ( !reference_ ) delete[] val;
}
};
std::ostream& operator << ( std::ostream& os, const IntegerMatrix& a );
class DoubleMatrix
{
private:
......@@ -436,10 +282,8 @@ class DoubleMatrix
void syevx(char uplo, std::valarray<double>& w, DoubleMatrix& z,
double abstol);
#if SCALAPACK
// permute the coeff of the matrix *this
void lapiv(char direc, char rowcol, IntegerMatrix& permutation);
#endif
void lapiv(char direc, char rowcol, const int *ipiv);
// compute eigenvalues (only) of symmetric matrix *this
// using the divide and conquer method of Tisseur and Dongarra
......@@ -698,10 +542,8 @@ class ComplexMatrix
// compute eigenvalues (only) of hermitian matrix *this
void heevd(char uplo, std::valarray<double>& w);
#if SCALAPACK
// permute the coeff of the matrix *this
void lapiv(char direc, char rowcol, IntegerMatrix &permutation);
#endif
void lapiv(char direc, char rowcol, const int *ipiv);
};
std::ostream& operator << ( std::ostream& os, const ComplexMatrix& a );
#endif
Markdown is supported
0% or .