使用MKL求解线性方程组的C++程序

#pragma once #include "StdAfx.h" #include <iostream> #include <math.h> #include "LinSolveLU.h" using namespace std; LinSolveLU::LinSolveLU(int N):n_(N),hcond_(0) { vec_=new LDOUBLE[n_]; w_=new MyMatrix<LDOUBLE>(n_,1); ipiv_=new int[n_]; iwork_=new int[n_]; work_=new LDOUBLE[4*n_]; } LinSolveLU::~LinSolveLU(void) { delete []ipiv_; delete []iwork_; delete []work_; delete []vec_; delete w_; } void LinSolveLU::Solve(const MyMatrix<LDOUBLE> &b,MyMatrix<LDOUBLE> &A)//利用LU分解求解非齐次线性方程组 { int info=0; int nrhs=1; LDOUBLE anorm=Norm1(A); w_->Copy(const_cast<MyMatrix<LDOUBLE> &>(b)); dgetrf(&n_,&n_,A.Data(),&n_,ipiv_,&info); dgecon("1",&n_,A.Data(),&n_,&anorm,&hcond_,work_,iwork_,&info); dgetrs("N",&n_,&nrhs,A.Data(),&n_,ipiv_,w_->Data(),&n_,&info); } LDOUBLE LinSolveLU::Norm1(const MyMatrix<LDOUBLE> &A) { LDOUBLE norm1=0; for(int i=0;i<n_;++i) { LDOUBLE temp=0; temp=cblas_dasum(n_,A.DataColumn(i+1),n_); vec_[i]=temp; } norm1=MaxVec(); return norm1; } LDOUBLE LinSolveLU::NormInf(const MyMatrix<LDOUBLE> &A) { LDOUBLE norm1=0; for(int i=0;i<n_;++i) { LDOUBLE temp=0; temp=cblas_dasum(n_,A.DataRow(i+1),1); vec_[i]=temp; } norm1=MaxVec(); return norm1; } LDOUBLE LinSolveLU::MaxVec() { LDOUBLE temp=0; for(int i=0;i<n_;++i) { if(vec_[i]>temp) temp=vec_[i]; } return temp; }