C++矩阵类(5)

double * pDbSrc=m.p; long nLen=m.width; int *is,*js,i,j,k; double d,p;

is = new int[nLen]; js = new int[nLen]; for(k=0;k

d=0.0;

for(i=k;i

p=fabs(pDbSrc[i*nLen + j]); if(p>d) {

d = p; is[k] = i; js[k] = j; } }

if(d+1.0==1.0) {

delete is; delete js; assert(0);

//return FALSE; }

if(is[k] != k)

for(j=0;j

p = pDbSrc[k*nLen + j];

pDbSrc[k*nLen + j] = pDbSrc[(is[k]*nLen) + j]; pDbSrc[(is[k])*nLen + j] = p; }

if(js[k] != k)

for(i=0; i

p = pDbSrc[i*nLen + k];

pDbSrc[i*nLen + k] = pDbSrc[i*nLen + (js[k])]; pDbSrc[i*nLen + (js[k])] = p; }

pDbSrc[k*nLen + k]=1.0/pDbSrc[k*nLen + k]; for(j=0; j

{

pDbSrc[k*nLen + j]*=pDbSrc[k*nLen + k]; }

for(i=0; i

for(j=0; j

pDbSrc[i*nLen + j] -= pDbSrc[i*nLen + k]*pDbSrc[k*nLen + j]; }

for(i=0; i

pDbSrc[i*nLen + k] *= -pDbSrc[k*nLen + k]; } }

for(k=nLen-1; k>=0; k--) {

if(js[k] != k)

for(j=0; j

p = pDbSrc[k*nLen + j];

pDbSrc[k*nLen + j] = pDbSrc[(js[k])*nLen + j]; pDbSrc[(js[k])*nLen + j] = p; }

if(is[k] != k)

for(i=0; i

p = pDbSrc[i*nLen + k];

pDbSrc[i*nLen + k] = pDbSrc[i*nLen +(is[k])]; pDbSrc[i*nLen + (is[k])] = p; } }

delete is; delete js; return m; }

Matrix Matrix::BiCubic(long n) {

assert(n>1); delete [] p; width=n; height=n;

p=new double[n*n];

//首行

*(p)=0.5; *(p+1)=2;

for(long k=2;k

for(long i=1;i

*(p+n*i+j)=0.5; j++;

*(p+n*i+j)=2; j++;

*(p+n*i+j)=0.5; j++; }

*(p+n*i+j)=0; } }

//末行

for(long l=0;l

return *this; }

Matrix Matrix::operator +(Matrix &m1) {

assert(m1.height==height && m1.width==width); long tmpHeight=m1.height; long tmpWidth=m1.width;

double * t=new double[tmpWidth*tmpHeight]; for(long i=0;i

*(t+tmpWidth*i+j)=*(m1.p+tmpWidth*i+j)+*(p+tmpWidth*i+j); } }

Matrix m(t,tmpHeight,tmpWidth); delete [] t; return m; }

Matrix Matrix::operator -(Matrix &m1) {

assert(m1.height==height && m1.width==width); long tmpHeight=m1.height; long tmpWidth=m1.width;

double * t=new double[tmpWidth*tmpHeight]; for(long i=0;i

*(t+tmpWidth*i+j)=*(p+tmpWidth*i+j)-*(m1.p+tmpWidth*i+j); } }

Matrix m(t,tmpHeight,tmpWidth); delete [] t; return m; }

Matrix Matrix::operator *(Matrix &m1) {

if(!this->isVector() && m1.isVector()){//左为数,右为矩阵 Matrix m; m=p[0]*m1; return m;

}else if(this->isVector() && !m1.isVector()){//左为矩阵,右为数 Matrix m;

m=*this*m1[0][0]; return m;

}else if(!this->isVector() && m1.isVector()){//左右都为数 double * t=new double[1]; t[0]=p[0]*m1[0][0]; Matrix m(t,1,1); delete [] t; return m;

}else if(this->isVector() && m1.isVector() && width==m1.height){//左为矩阵,右为矩阵 double sum;

double * t=new double[height*m1.width]; for(long i=0;i

for(long k=0;k

sum+=(*(p+width*i+k))*(m1[k][j]); }

*(t+m1.width*i+j)=sum; } }

Matrix m(t,height,m1.width); delete [] t; return m; }else{

assert(0);//未知运算 return *this;

} }

Matrix operator*(double alpha,Matrix & m1) {

Matrix m=m1;

for(long i=0;i

return m; }

Matrix Matrix::operator*(double alpha) {

return alpha*(*this); }

Matrix Matrix::operator+=(Matrix & m) {

return *this+m; }

Matrix Matrix::operator-=(Matrix & m) {

return *this-m; }

Matrix Matrix::operator *=(double alpha) {

return *this*alpha; }

Matrix Matrix::operator *=(Matrix & m1) {

return *this*m1; }

Matrix sqrt(Matrix m) {

assert(m[0][0]>=0 && m.getHeight()==1 && m.getWidth()==1); m[0][0]=sqrt(m[0][0]); return m; }

double abs(Matrix & m) {

double sum=0;

for(long i=0;i