www.pudn.com > roadextr.rar > Matrix.cpp

#include	"stdafx.h" 
#include	"valdef.h" 
#include	 
#include	 
#include	"Matrix.h" 
 
void MatrixMultiply(double *mA, double *mB, int rowa, int cola, int colb, double *mC) 
{ 
	int	i, j, k; 
 
	for (i=0 ; i 2)  
  { /*if the matrix is 3 x 3 or bigger */ 
	for (i=0;it2) 
//				if (fabs(a[i-1][j-1])>fabs(w)) 
				{ 
					w = *(a + (i-1)*n + j-1);// w = a[i-1][j-1]; 
					i0 = i; 
					j0 = j; 
				} 
			} 
		t2 = w; 
		if (w<0) t2 = -w; 
		if (false) /// t2<=EP 
//		if (fabs(w)<=EP) 
		{ 
			delete []p; 
			delete []q; 
			delete []b; 
			delete []c; 
			return false; 
		} 
		if (i0!=k) 
			for (j=1;j<=n;j++) 
			{ 
				z = *(a + (i0-1)*n + j-1);/// z = a[i0-1][j-1]; 
				*(a + (i0-1)*n + j-1) = *(a + (k-1)*n + j-1); /// a[i0-1][j-1] = a[k-1][j-1]; 
				*(a + (k-1)*n + j-1) = z;/// a[k-1][j-1] = z; 
			} 
		if (j0!=k) 
			for (i=1;i<=n;i++) 
			{ 
				z = *(a + (i-1)*n + j0-1);/// z = a[i-1][j0-1]; 
				*(a + (i-1)*n + j0-1) = *(a + (i-1)*n + k-1); /// a[i-1][j0-1] = a[i-1][k-1]; 
				*(a + (i-1)*n + k-1) = z;//a[i-1][k-1] = z; 
			} 
		p[k-1] = i0; 
		q[k-1] = j0; 
		for (j=1;j<=n;j++) 
		{ 
			if (j==k) 
			{ 
				b[j-1] = 1/w; 
				c[j-1] = 1; 
			} 
			else 
			{ 
				b[j-1] = -*(a + (k-1)*n + j-1)/w;//// -a[k-1][j-1]/w; 
				c[j-1] = *(a + (j-1)*n + k-1);//// a[j-1][k-1]; 
			} 
			*(a + (k-1)*n + j-1) = *(a + (j-1)*n + k-1) = 0.0;/// a[k-1][j-1] = a[j-1][k-1] = 0.0; 
		} 
		for (i=1;i<=n;i++) 
			for (j=1;j<=n;j++) 
				/// a[i-1][j-1] = a[i-1][j-1] + c[i-1]*b[j-1]; 
				*(a + (i-1)*n + j-1) = *(a + (i-1)*n + j-1) + c[i-1]*b[j-1]; 
	} 
    for (k=n;k>=1;k--) 
    { 
        i0 = p[k-1]; 
        j0 = q[k-1]; 
        if (i0!=k) 
			for (i=1;i<=n;i++) 
			{ 
				z = *(a + (i-1)*n + i0-1);/// a[i-1][i0-1]; 
				*(a + (i-1)*n + i0-1) = *(a + (i-1)*n + k-1);//a[i-1][i0-1] = a[i-1][k-1]; 
				*(a + (i-1)*n + k-1) = z;/// a[i-1][k-1] = z; 
			} 
		if (j0!=k) 
			for (j=1;j<=n;j++) 
			{ 
				z = *(a + (j0-1)*n + j-1);/// a[j0-1][j-1]; 
				*(a + (j0-1)*n + j-1) = *(a + (k-1)*n + j-1);/// a[j0-1][j-1] = a[k-1][j-1]; 
				*(a + (k-1)*n + j-1) = z; // a[k-1][j-1] = z; 
			} 
	} 
	delete[]p; 
	delete[]q; 
    delete[]b; 
    delete[]c; 
    return true; 
}