I meet a serious problem when I was coding my programme. Firstly, I will state my implementation.
I have a Conjugate Gradient Implementation to solve linear system. In my densematrixvector.h, I have two classes, denseMatrix and denseVector, and I put these two arrays into the function template < class M, class V > int ConjugateGradient(const M &A, const V &b, V &x, const int &maxstep, const int &eps), here, M is denseMatrix, V is denseVector.
In function ConjugateGradient, I am using a function void gemv(const double &alpha, const DenseMatrix &A, const DenseVector &x, const double &beta, DenseVector &y) to solve the linear system.
But this code did not take advantage of the sparsity of the matrix. So I want to modify the code so that it does. In particular, the matrix type passed to the ConjugateGradient function represent dense matrix. I want to modify the code to pass it sparse matrix in ccs format. this class transfers information of sparse matrix. So I modify the declaration of function gemv like this:
template <class M, class V >void gemv(const double &alpha, const M &A, const V &x, const double &beta, V &y)
And the declaration of matCCS is template < class T >class matCCS.
finally, my computer tells me that undefined reference to.
and the code I attach below.
Please help me. Thank you in advance.
mainconjugategradientmethodmatCCS.cc
#include "densematrixvector.h"
#include "ConjugateGradient.h"
#include <iostream>
#include "sparseMatrix.h"
int main(){
matCOO<double> Acoo;
readMatrixMarket("small.mtx", Acoo);
matCCS<double> A; // constuct a CCS matrix
COOtoCCS(Acoo, A); // convert the COO matrix Acoo to CCS A
int m, n;
m=A.m;
n=A.n;
DenseVector x;
DenseVector b;
int nx = 1;
double * xt = new double [m*nx];
for (int i = 0; i < m ;++i) xt[i] =1.;
x.a= xt;
x.n = n;
int nrhs = 1;
double * bt = new double [m*nrhs];
for (int i = 0; i < m ;++i) bt[i] =0.;
b.a = bt;
b.n = n;
int nstep = ConjugateGradient(A, b, x, 100, 0.001);
if (nstep > 0){
std::cout << "converged in nstep " <<nstep << std::endl;
std::cout<< " solution is " << x.a[0] << " " << x.a[1] << std::endl;
}
else {
std::cout << "CG could not converge " <<nstep << std::endl;
}
}
densematrixvector.cc
#include "densematrixvector.h"
#include <algorithm>
#include <iostream>
DenseVector::DenseVector():a(0), n(0), owner(false){};
DenseVector::DenseVector(const int & _n):a(new double[_n]), n(_n), owner(true)
{
std::fill(a,a+n, 0.);
};
DenseVector::DenseVector(const DenseVector &in):a(new double[in.n] ), n(in.n), owner(true){
std::copy(in.a, in.a+n, a);
};
DenseVector::~DenseVector(){
if(owner) delete []a;
owner =false;
};
int size(const DenseVector &in){
return in.n;
}
void copy(const DenseVector &x, DenseVector &y){
const int n = size(x);
if (n!=size(y)){ std::cout << " Size error in copy " << __FILE__ <<
__LINE__ << std::endl;}
std::copy(x.a, x.a+n, y.a);
}
double dot(const DenseVector & x, const DenseVector & y){
const int n = size(x);
if (n!=size(y)){ std::cout << " Size error in dot " << __FILE__ << __LINE__ << std::endl;}
return ddot_(n, x.a, 1, y.a, 1);
}
void scal(const double& alpha, DenseVector &x){
dscal_(size(x), alpha, x.a, 1);
}
void axpy( const double &alpha, const DenseVector & x, DenseVector &y){
const int n = size(x);
if (n!=size(y)){ std::cout << " Size error in axpy " << __FILE__ <<
__LINE__ << std::endl; throw;}
daxpy_(n, alpha, x.a, 1, y.a, 1 );
}
template <class M, class V>
void gemv(const double& alpha, const M &A, const V &x, const double &beta, V &y )
int m =A.m;
int n =A.n;
int nnz=A.nnz;
if ((size(x) != n) || (size(y) != m)) { std::cout << " Size error in gemv "` `<< __FILE__ << __LINE__ << std::endl; throw;}
//dgemv_('N', m,n, alpha, A.a, m, x.a, 1, beta, y.a, 1 );
for(int i=0; i<m; i++) {
y.a[i]=A.a[i]*x.a[A.lineindex[i]];
}
}
densematrixvector.h
#ifndef _densematrixvector_
#define _densematrixvector_
...
int size(const DenseVector &in);
void copy(const DenseVector &x, DenseVector &y );
double dot(const DenseVector & x, const DenseVector & y);
void scal(const double& alpha, DenseVector &x);
void axpy( const double &alpha, const DenseVector & x, DenseVector &y);
//void gemv(const double& alpha, const DenseMatrix &A, const DenseVector &x,
const double &beta, DenseVector & y );
template <class M, class V>
void gemv(const double& alpha, const M &A, const V &x, const double &beta, V &y );
...
#endif
ConjugateGradient.h
#include <iostream>
template <class M, class V>
int ConjugateGradient(const M & A, const V& b, V &x, const int & maxstep, const int &eps)
{
int n = size(b);
std::cout << n << std::endl;
//r = b - A x;
V r(b);
gemv(-1., A, x, 1., r);
...
The report error is attached below
[master@localhost part3ddis]$ g++ maintestconjugategradientmatCCS.cc densematrixvector.cc mmio.cc libsuperlu_4.3.a libblas.a -lgfortran
/tmp/ccg7iPUj.o: In function `int ConjugateGradient<DenseMatrix, DenseVector>(DenseMatrix const&, DenseVector const&, DenseVector&, int const&, int const&)':
maintestconjugategradientmatCCS.cc:(.text._Z17ConjugateGradientI11DenseMatrix11DenseVectorEiRKT_RKT0_RS5_RKiSA_[_Z17ConjugateGradientI11DenseMatrix11DenseVectorEiRKT_RKT0_RS5_RKiSA_]+0xb2): undefined reference to `void gemv<DenseMatrix, DenseVector>(double const&, DenseMatrix const&, DenseVector const&, double const&, DenseVector&)'
maintestconjugategradientmatCCS.cc:(.text._Z17ConjugateGradientI11DenseMatrix11DenseVectorEiRKT_RKT0_RS5_RKiSA_[_Z17ConjugateGradientI11DenseMatrix11DenseVectorEiRKT_RKT0_RS5_RKiSA_]+0x208): undefined reference to `void gemv<DenseMatrix, DenseVector>(double const&, DenseMatrix const&, DenseVector const&, double const&, DenseVector&)'
collect2: error: ld returned 1 exit status
