THE Matrix

// Copyright (C) Liang Yi, 2004

// Finally modified: 2004.8. Hangzhou

//Vector3d.h

#ifndef Vector3D_Header

#define vector3d_header

#include

#include

Using namespace std;

Template

Class vector3d

{

PUBLIC:

Vector3D (t _x = 0, t _y = 0, t_z = 0.0)

{

X = _x;

Y = _y;

z = _z;

}

Vector3D (Const Vector3D & V)

{

x = v.x;

y = v.y;

z = v.z;

}

Template

VECTOR3D (Const in PV)

{

X = PV [0];

Y = PV [1];

Z = PV [2];

}

Inline T Magnitude ()

{

Return SQRT (x * x y * y z * z);

}

Inline void normalize ()

{

T MAGNITUDE = MAGNITUDE ();

IF (! magnetude) return;

X / = magnise;

Y / = magnise;

z / = magnise;

}

Inline Vector3D getnormalizedVector ()

{

T MAGNITUDE = MAGNITUDE ();

IF (! Magnitude) Return * this;

Return (* this) / magnitude;

}

Inline Vector3D Operator (Vector3D & V)

{

Return VECTOR3D (x v.x, y v.y, z v.z);

}

Inline Vector3D Operator - (Vector3D & V)

{

Return VECTOR3D (x-v.x, y-v.y, z-v.z);

}

Inline T Operator * (Vector3D & V)

{

Return X * v.x y * v.y z * v.z;

}

Inline Vector3D x (Vector3D & V)

{

Return VECTOR3D (Y * V.x-x * v.y, z * v.x-x * v.z, x * v.y-y * z.x);

}

Inline Vector3D Operator (t t)

{

Return VECTOR3D (x t, y t, z t);

}

Inline Vector3D Operator - (t t)

{

Return VECTOR3D (x-T, Y-T, Z-T);

}

Inline Vector3D Operator * (t t)

{

Return VECTOR3D (x * t, y * t, z * t);

}

Inline Vector3D Operator / (t T)

{

Return VECTOR3D (x / t, y / t, z / t);

}

T & Operator [] (size_t index) // Nonconst, But Without & (Can Be Left-Value) {

Switch (INDEX)

{

Case 0:

Return X;

Case 1:

Return Y;

DEFAULT:

Return Z;

}

}

Friend Ostream & Operator << (Ostream & OS, Const Vector3D & V)

{

Return OS << v.x << "/ t" << v.y << "/ t" << v.z << Endl;

}

T x, y, z;

}

#ENDIF

//matrix.h

#ifndef matrix_h

#define matrix_h

#include "vector3d.h"

#include

#include

Using namespace std;

Template

Class Matrix

{

PUBLIC:

Template

Matrix (const in pdata, size_t nr, size_t nc): NROW (NR), NCOL (NC) // PDATA continuous container (array)

{

T * p = (t *) (& (PDATA [0]));

Copy (p, p nr * nc, insert))); // data.reserve (NR * NC); COPY (P, P NR * NC, DATA.BEGIN ()); Faster

}

Matrix ()

{

Ncol = 0; nROW = 0;

}

Matrix (size_t _nrow, size_t _ncol = 0)

{

IF (_ncol == 0) _ncol = _nrow;

nrow = _nrow;

Ncol = _ncol;

IF (NROW> 0)

Identity ();

}

void identity ()

{

Data.resize (nROW * NCOL, 0);

For (size_t i = 0; i <(ncol? nrow: ncol); i )

{

Data [i * ncol i] = 1.0;

}

}

Friend Ostream & Operator << (Ostream & OS, Const Matrix & M)

{

IF (m.nrow <= 0 || m.ncol <= 0) Return OS;

For (SIZE_T I = 0; I

{

For (SIZE_T J = 0; J

{

OS << M.Data [i * m.ncol j] << "/ t";

}

OS << Endl;

}

Return OS;

}

Void SetValue (t v, size_t r, size_t c)

{

Data [t * nrow c] = v;

}

T getValue (size_t r, size_t c) const

{

Return Data [T * NROW C];

}

Void getRotor (Matrix & OUT)

{

Out.Data.Rserve (data.size ());

Out.ncol = nrow;

Out.nrow = ncol; for (size_t i = 0; i

{

For (SIZE_T J = 0; J

{

Out.data [i * nrow j] = data [j * ncol i];

}

}

}

Matrix getRotor ()

{

Matrix OUT;

Out.data.resize (data.size (), 0);

Out.ncol = nrow;

Out.nrow = ncol;

For (size_t i = 0; i

{

For (SIZE_T J = 0; J

{

Out.data [i * nrow j] = data [j * ncol i];

}

}

Return Out;

}

Void ExchangeRowPOS (SIZE_T R1, SIZE_T R2)

{

IF (r1> = nrow || r2> = nrow || r2 == r1) Return;

T TMP;

For (size_t i = 0; i

{

TMP = DATA [R1 * NCOL I];

Data [R1 * NCOL I] = DATA [R2 * NCOL I];

DATA [R2 * NCOL I] = TMP;

}

}

Void ExchangeColpos (SIZE_T C1, SIZE_T C2)

{

IF (c1> = nrow || c2> = nrow || c2 == c1) return;

T TMP;

For (SIZE_T I = 0; i

{

TMP = DATA [i * nrow c1];

Data [i * nrow c1] = data [i * nrow c2];

Data [i * nrow c2] = TMP;

}

}

Vector getX (Vector b)

{

Vector x;

x.resize (nrow, 0);

IF (NROW <2 || ncol <2 || ncol! = nrow || b.size ()! = nrow) Return X;

Matrix m = * this;

SIZE_T I, J, K;

FOR (k = 0; k

{

IF (M.Data [K * Ncol K] == 0)

{

FOR (i = k 1; i

{

IF (M.Data [i * ncol k]! = 0)

{

M.ExchangerowPos (K, I);

B [K] = B [I];

B [i] = b [k] -b [i];

B [k] = b [k] -b [i];

Break;

}

}

IF (i == nrow)

Return X;

}

FOR (i = k 1; i

{

T rate = m.data [i * ncol k] /m.data [k * ncol k];

For (j = k; j

{

M.DATA [i * ncol j] - = rate * m.data [k * ncol j];

}

B [I] - = rate * b [k];

}

}

X [nrow-1] = b [nrow-1] /m.data [m.data.size () - 1]; for (int A = (int) nrow-2; a> = 0; A -)

{

T tmpsum = 0;

For (int b = (int) ncol-1; b> a; b -)

{

Tmpsum = x [b] * (M.Data [a * ncol b]);

}

X [A] = (B [a] -tmpsum) /m.data [a * ncol a];

}

Return X;

}

T getModel ()

{

IF (NROW <2 || ncol <2 || ncol! = nrow) return 0;

Matrix m = * this;

SIZE_T I, J, K;

FOR (k = 0; k

{

IF (M.Data [K * Ncol K] == 0)

{

FOR (i = k 1; i

{

IF (M.Data [i * ncol k]! = 0)

{

M.ExchangerowPos (K, I);

Break;

}

}

IF (i == nrow)

Return 0;

}

FOR (i = k 1; i

{

T rate = m.data [i * ncol k] /m.data [k * ncol k];

For (j = k; j

{

M.DATA [i * ncol j] - = rate * m.data [k * ncol j];

}

}

}

T TMP = 1;

For (k = 0; k

{

TMP * = M.Data [k * ncol k];

}

Return TMP;

}

Void getNegativeMatrix (Matrix & Out)

{

IF (NROW <2 || NCOL <2 || ncol! = nrow) return;

Out.data.clear ();

OUT = * this;

SIZE_T K, I, J;

Vector rc_list;

Size_T RCI, RCJ;

T TMP;

FOR (k = 0; k

{

// Select the main yuan, stable, can remove this part with stability

TMP = FABS (Out.Data [k * ncol k]);

RCI = K;

RCJ = K;

For (i = k; i

{

For (j = k; j

{

IF (Fabs (Out.Data [i * ncol j])> TMP)

{

RCI = I;

RCJ = J;

TMP = FABS (Out.Data [i * ncol j]);

}

}

}

/ / Save the ranks exchange information when the selection

IF (RCI! = K)

{

Out.exchangerowpos (K, RCI);

Row_col rc;

Rc.RC1 = K;

Rc.RC2 = RCI;

rc.brow = true;

RC_List.push_back (rc);

}

IF (RCJ! = K)

{

Out.exchangecolpos (k, rcj);

Row_col rc;

Rc.RC1 = K;

Rc.RC2 = RCJ;

Rc.brow = false;

RC_List.push_back (rc);

}

// select the main element end

Out.data [k * ncol k] = 1 / OUT.DATA [K * NCOL K];

For (j = 0; j

{

IF (k! = j)

Out.data [k * ncol j] = out.data [k * ncol j] * out.data [k * ncol k];

}

For (i = 0; i

{

For (j = 0; j

{

IF (k! = j && k! = i)

Out.data [i * ncol j] = out.data [i * ncol j] -Out.data [i * ncol k] * out.data [k * ncol j];

}

}

For (i = 0; i

{

IF (k! = i)

Out.data [i * ncol k] * = - out.data [k * ncol k];

}

}

/ / Restore the ranks exchange

For (k = rc_list.size (); k> 0; k -)

{

i = k-1;

IF (rc_list [i] .brow)

{

Out.exchangecolpos (rc_list [i] .rc1, rc_list [i] .rc2);

}

Else

{

OUT.EXCHANGEROWPOS (RC_List [i] .rc1, rc_list [i] .rc2);

}

}

}

Void getJoint (Matrix & Out)

{

Out.data.clear ();

GetNegativeMatrix (OUT);

Vector :: item ney = out.data.begin ();

T m = getModel ();

For (; it! = out.data.end (); iTer )

{

(* iter) * = m;

}

}

Matrix Operator (Const Matrix & M)

{

Matrix OUT = * this;

For (SIZE_T I = 0; i

{

Out.data [i] = data [i] m.data [i];

}

Return Out;

}

Matrix Operator - (Const Matrix & M)

{

Matrix OUT = * this;

For (SIZE_T I = 0; i

{

Out.data [i] = data [i] -m.data [i];

}

Return Out;

}

Matrix Operator * (Const Matrix & M)

{

IF (ncol! = m.nrow) Return Matrix ();

Vector OUT;

Out.resize (nrow * m.ncol, 0);

SIZE_T I, J, K;

For (i = 0; i

{

For (j = 0; j

{

For (k = 0; k

OUT [i * m.ncol j] = data [i * ncol k] * m.data [k * m.ncol j];

}

Return Matrix (OUT, NROW, M.NCOL);

}

Void Add (Const Matrix & M)

{

For (SIZE_T I = 0; i

{

Data [i] = m.data [i];

}

}

Void Sub (Const Matrix & M)

{

For (SIZE_T I = 0; i

{

Data [i] - = m.data [i];

}

}

Void Multiply (Const Matrix & M)

{

IF (ncol! = m.nrow) return;

Vector OUT;

Out.resize (nrow * m.ncol, 0);

SIZE_T I, J, K;

For (i = 0; i

{

For (j = 0; j

{

For (k = 0; k

OUT [i * m.ncol j] = data [i * ncol k] * m.data [k * m.ncol j];

}

}

Data = OUT;

Ncol = m.ncol;

}

Vector3D Multiply (Const Vector3D & V, T R = 1.0)

{

IF (ncow! = 4 || ncol! = 4) Return Vector3D ();

SIZE_T I, K;

T out [4] = {0,0,0,0};

T in [4] = {v.x, v.y, v.z, r};

For (i = 0; i <4; i )

{

For (k = 0; k <4; k )

OUT [i] = DATA [i * ncol k] * in [k];

}

Return Vector3D (OUT);

}

///

Vector DATA;

SIZE_T NROW, NCOL;

Private:

Struct row_col

{

SIZE_T RC1;

SIZE_T RC2;

Bool brow;

}

}

#ENDIF