Linux下采用lapack科学计算包来实现二维矩阵

1、调用lapack包的方式，写一个头文件 my_lapack.h （Linux如何安装lapack包？）

#ifndef MY_LAPACK_H
#define MY_LAPACK_H

extern "C"
{
	#include "lapacke.h"
	extern void dgesv_(int*,int*,double*,int*,int*,double*,int*,int*);
	extern void sgesv_(int*,int*,float*,int*,int*,float*,int*,int*);
}

#endif //MY_LAPACK_H

2、使用C++来实现二维矩阵，并采用模板类。具体如下：my_matrix.h

#ifndef MY_MATRIX_H
#define MY_MATRIX_H

#include <iostream>
#include "my_lapack.h"
#include <string>

template <typename T>
class my_matrix{

  //my_matrix类的成员变量
  int num_rows;
  int num_columns;
  T** ptr;
  std::string name;

  public:
  //my_matrix类的默认构造函数
  my_matrix();

  // 自定义的构造函数
  my_matrix(int Nrows, int Ncols, std::string mat_name);

  // 拷贝构造函数
  my_matrix(const my_matrix<T>& mat);

  // my_matrix类的公共成员函数display()
  void display();

  // my_matrix类的析构函数
  ~my_matrix();

  // =运算符重载函数
  my_matrix<T>& operator=(const my_matrix<T>& mat);

  //()运算符重载函数
  T& operator()(int row, int col);

  //+运算符重载函数
  my_matrix<T> operator+(const my_matrix<T>& mat);

  //|运算符重载函数
  my_matrix<T> operator|(const my_matrix<T>& B);

};

//my_matrix类的默认构造函数
template <typename T>
my_matrix<T>::my_matrix(){
  num_rows = 2;
  num_columns = 2;
  name = "default";
  //ptr变量的初始化操作
  ptr = new T*[num_columns];
  ptr[0] = new T[num_columns * num_rows];
  for (int jj = 1; jj < num_columns; jj++){
    ptr[jj] = ptr[0] + num_rows * jj;
  }
}

// user constructor
template <typename T>
my_matrix<T>::my_matrix(int Nrows, int Ncols, std::string mat_name){
  try{
    if (Nrows < 1 || Ncols < 1){
      throw "matrix dimensions do not make sense.";
    }
    else {
        num_rows = Nrows;
        num_columns = Ncols;
        name = mat_name;
        ptr = new T*[num_columns];
        ptr[0] = new T[num_columns * num_rows];
        for (int jj = 1; jj < num_columns; jj++){
            ptr[jj] = ptr[0] + num_rows * jj;
        }
    }
  }
  catch (const char* problem){
    std::cout << "exception encountered: " << problem << std::endl;
    ptr = new T*[1];
    ptr[0] = new T[1];
  }
}

// copy constructor
template <typename T>
my_matrix<T>::my_matrix(const my_matrix<T>& mat){
  num_rows = mat.num_rows;
  num_columns = mat.num_columns;
  name = mat.name;
  // initialize with column major format
  ptr = new T*[num_columns];
  ptr[0] = new T[num_columns * num_rows];
  for (int ii = 1; ii < num_columns; ii++){
    ptr[ii] = ptr[0] + num_rows * ii;
  }
  for (int jj = 0; jj < num_columns; jj++){
    for (int ii = 0; ii < num_rows; ii++){
      ptr[jj][ii] = mat.ptr[jj][ii];
    }
  }
}

// destructor
template <typename T>
my_matrix<T>::~my_matrix(){
  delete[] ptr[0];
  delete[] ptr;
}

// for printing out elements
template <typename T>
void my_matrix<T>::display(){
  std::cout << this->name << " = \n";
  for (int ii = 0; ii < this->num_rows; ii++){
    for (int jj = 0; jj < this->num_columns; jj++){
      std::cout << this->ptr[jj][ii] << " ";
      //std::cout << (*this)(ii,jj) << " ";
    }
    std::cout << "\n";
  }
  std::cout << "\n";
}

// implementation of operators
template <typename T>
my_matrix<T>& my_matrix<T>::operator=(const my_matrix<T>& mat){
  if (this != &mat){
    // deallocate memory
    for (int jj = 0; jj < num_columns; jj++){
      delete[] ptr[jj];
    }
    delete[] ptr;

    // reallocate memory
    this->num_rows = mat.num_rows;
    this->num_columns = mat.num_columns;
    this->name = mat.name;
    // initialize with column major format
    ptr = new T*[this->num_columns];
    ptr[0] = new T[this->num_columns * this->num_rows];
    for (int jj = 1; jj < this->num_columns; jj++){
        ptr[jj] = ptr[0] + this->num_rows * jj;
    }
    // copy elements
    for (int jj = 0; jj < this->num_columns; jj++){
        for (int ii = 0; ii < this->num_rows; ii++){
            ptr[jj][ii] = mat.ptr[jj][ii];
        }
    }
  }

  return *this;
}

//该()运算符重载按照坐标来直接进行获取值，如(0,1)则我们直接使用(0,1)来访问这个位置
template <typename T>
T& my_matrix<T>::operator()(int row, int col)
{
    return ptr[row][col];
}

template <typename T>
my_matrix<T> my_matrix<T>::operator+(const my_matrix<T>& mat){
  my_matrix<T> sum(this->num_rows, this->num_columns, "sum");
   for (int jj = 0; jj < this->num_columns; jj++){
        for (int ii = 0; ii < this->num_rows; ii++){
            sum.ptr[jj][ii] = ptr[jj][ii] + mat.ptr[jj][ii];
        }
    }
  return sum;
}

//|来实现二维矩阵的逆矩阵！ 采用的是科学计算包lapack来实现的！~
template <>
my_matrix<double> my_matrix<double>::operator|(const my_matrix<double>& B){

  int N = this->num_rows;
  int NRHS = B.num_columns;
  my_matrix<double> A_copy(*this);
  my_matrix<double> B_copy(B);
  int LDA = this->num_rows;
  int* IPIV = new int[N];
  int LDB = B.num_rows;
  int INFO = 0;

  dgesv_(&N, &NRHS, A_copy.ptr[0], &LDA, IPIV, B_copy.ptr[0], &LDB, &INFO);

  std::cout << "INFO = " << INFO << std::endl;

  delete[] IPIV;
  return B_copy;
}

template <>
my_matrix<float> my_matrix<float>::operator|(const my_matrix<float>& B){

  //continue here for sgesv_
  int N = this->num_rows;
  int NRHS = B.num_columns;
  my_matrix<float> A_copy(*this);
  my_matrix<float> B_copy(B);
  int LDA = this->num_rows;
  int* IPIV = new int[N];
  int LDB = B.num_rows;
  int INFO = 0;

  sgesv_(&N, &NRHS, A_copy.ptr[0], &LDA, IPIV, B_copy.ptr[0], &LDB, &INFO);

  std::cout << "INFO = " << INFO << std::endl;

  delete[] IPIV;
  return B_copy;
}

#endif

3、主函数里如何调用这些类函数，实现二维矩阵！main.cpp

#include "my_matrix.h"
#include <iostream>
#include "my_lapack.h"
#include <string>

//测试这些接口
void test_Operators_function(void)
{
	my_matrix<double> testA(2, 2, "A");
	my_matrix<double> testB(2, 2, "B");
	testA(0,0) = 1;
    testA(0,1) = 1;
    testA(1,0) = 1;
    testA(1,1) = 1;

    testB(0,0) = 2;
    testB(0,1) = 2;
	testB(1,0) = 2;
    testB(1,1) = 2;
	my_matrix<double> testsum = testA + testB;
	testsum.display();
}
int main(void){
	test_Operators_function();
  return 0;
}

4、我们使用make工具来管理这些文件，Makefile文件内容如下：

cc=g++
obj=main.o
target=main
cflags=-llapacke -llapack -lrefblas -lgfortran -lm -w

$(target):$(obj)
	$(cc)  $^  -o  $@  $(cflags)

main.o:main.cpp lapack_prototypes.h my_matrix.h
	$(cc)   -c  $<   -o  $@ 

clean:
	rm $(obj)  $(target)