Returns the inverse of the square matrix.

I wrote a function to calculate the inverse of the square matrix, because it is not a pseudo inverse of the square matrix, and the function should be written again when the append code is available.

I need matrixdet. h. I don't know how to find it in my csdnblog. I have not classified it, so I need to find the function code or something on my csdnblog.

/*************************************** ************************************
Copyright (c), all, begtostudy. NERC. zzu.
File Name: matrix_invers.h
Author: begtostudy version: Date: 2006.11.20
Description:
Others:
Function list:
History:
1. Date:
Author:
Modification:
**************************************** ***********************************/

// Put it in the header file matrix_invers.h. It will be used in my articles in the future. Don't say you don't know.
# Include <math. h>
# Include <memory. h>
# Include <vector>
Namespace jks
{
//////////////////////////////////////// //////////////////////////////////

Template <typename T>
Double * col1_( T * a, T * B, int N, double * & X)
{
Int I, j, T, K;
Double P;
If (x = NULL)
{
X = new double [N];
}
# Ifdef _ WINBASE _
Else
Assert (! : Isbadwriteptr (x, N * sizeof (double )));
Assert (! : Isbadreadptr (A, N * n * sizeof (t )));
Assert (! : Isbadreadptr (B, n * sizeof (t )));
# Endif

T * c = new T [N * (n + 1)];
For (I = 0; I <n; I ++)
{
Memcpy (C + I * (n + 1), a + I * n, n * sizeof (t ));
* (C + I * (n + 1) + n) = * (B + I );//!
}

For (I = 0; I <= n-2; I ++)
{
K = I;
For (j = I + 1; j <= n-1; j ++)
If (FABS (* (C + J * (n + 1) + I)> (FABS (* (C + K * (n + 1) + I ))))
K = J;
If (K! = I)
For (j = I; j <= N; j ++)
{
P = * (C + I * (n + 1) + J );
* (C + I * (n + 1) + J) = * (C + K * (n + 1) + J );
* (C + K * (n + 1) + J) = P;
}
For (j = I + 1; j <= n-1; j ++)
{
P = (* (C + J * (n + 1) + I)/(double) (* (C + I * (n + 1) + I ));
For (t = I; t <= N; t ++)
* (C + J * (n + 1) + T)-= p * (C + I * (n + 1) + T ));
}
}
For (I = n-1; I> = 0; I --)
{
For (j = n-1; j> = I + 1; j --)
(* (C + I * (n + 1) + n)-= x [J] * (C + I * (n + 1) + j ));
X [I] = * (C + I * (n + 1) + N)/(double) (* (C + I * (n + 1) + I ));
}

Delete [] C;
 
Return X;
}
Template <class T>
Void verticalcpy (T * Top, int towidth, const T * fromp, int fromheight, int fromwidth = 1)
{
// Assert (top! = NULL );
// Assert (fromp! = NULL );
// Assert (fromheight> 0 );

If (fromwidth = towidth)
{
STD: Copy (fromp, fromp + fromwidth * fromheight, top );
Return;
}
 
Int I, J;
For (I = 0; I <fromheight; I ++)
{
T * PT = Top + I * towidth;
For (j = 0; j <fromwidth; j ++)
{
* PT ++ = * fromp ++;
}
}
}
Template <typename T>
T * matix_reversal (T * P, int N)
{
// Assert (! : Isbadreadptr (p, N * n * sizeof (t )));

T * PTO = new T [N * n];
Int I;
For (I = 0; I <n; I ++)
{
Verticalcpy (PTO + I, N, P + I * n, n );
}

Memcpy (p, PTO, N * n * sizeof (t ));
Return P;

}
Template <class T>
T * matrix_multiplication (T * LPA, T * LPB, T * LPC,
Unsigned short M, unsigned short S, unsigned short N)
{
// Assert (LPC! = LPA & LPC! = LPB );
// Assert (LPC! = NULL );
// Assert (M! = 0 & n! = 0 & S! = 0 );
// Assert (! : Isbadwriteptr (LPC, M * n * sizeof (t )));
// Assert (! : Isbadreadptr (LPA, M * S * sizeof (t )));
// Assert (! : Isbadreadptr (LPB, S * n * sizeof (t )));

Unsigned short I, J, K;
// For (I = 0; I <m; I ++)
// For (j = 0; j <n; j ++)
// * (LPC + I * n + J) = 0;
Memset (LPC, 0, M * n * sizeof (t ));

For (I = 0; I <m; I ++)
{
For (j = 0; j <n; j ++)
{
For (k = 0; k <s; k ++)
{
* (LPC + I * n + J) + = * (LPA + I * s + k) ** (LPB + K * n + J );
// Cerr <* (LPA + I * s + k) <'/t' <* (LPB + K * n + J) <Endl;
}
}
}

Return LPC;
}
Inline double * matrix_inverse (double * P, int N, double * PTO)
{
// Assert (! : Isbadreadptr (p, N * n * sizeof (double )));

If (FABS (jks: matrix_det (p, n) <1e-20)
{
Return NULL;
}

Int I;

If (PTO = NULL)
{
PTO = new double [N * n];
}
// Else
// Assert (! : Isbadwriteptr (PTO, N * n * sizeof (double )));

Double * E = new double [N];

For (I = 0; I <n; I ++)
{
Memset (E, 0, N * sizeof (double ));
E [I] = 1;
Double * t = PTO + I * N;
Col1_( P, E, N, t );
}
Delete [] E;

Matix_reversal (p, N );

Return PTO;
}
//////////////////////////////////////// //////////////////////////////////
}

========================== Example

Void test ()
{
Double PP [3] [3] = {1, 2, 1, 1, 2, 2, 1, 4, 2 };
Double * P = PP [0];
// Jks: dump (p, 3, 3 );

Double * TP = new double [3*3];
If (jks: matrix_inverse (p, 3, TP) = NULL)
{
Cerr <"erro! /N ";
Return;
}

// Jks: dump (TP, 3, 3 );

Double * Rp = new double [3*3];
Jks: matrix_multiplication (p, TP, RP, 3, 3 );

// Jks: dump (RP, 3, 3 );

Int I, J;
For (I = 0; I <3; I ++)
{
For (j = 0; j <3; j ++)
{
Cerr <RP [I * 3 + J] <'';
} Cerr <Endl;
}

Delete [] RP;
}