Đến nội dung

sans_amour

sans_amour

Đăng ký: 13-11-2006
Offline Đăng nhập: 11-03-2008 - 00:41
-----

Trong chủ đề: GraphCalc 4.0.1

06-10-2007 - 15:07

Merci Magus ;:D

Trong chủ đề: Nghịch đảo ma trận bằng các ngôn ngữ

24-08-2007 - 01:00

lu.bas:

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'This code is generated by the AlgoPascal translator
'
'This code is distributed under the ALGLIB license
'	(see http://www.alglib.net/copyrules.php for details)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'LU decomposition of a general matrix of size MxN
'
'The subroutine calculates the LU decomposition of a rectangular general
'matrix with partial pivoting (with row permutations).
'
'Input parameters:
'	A   -   matrix A whose indexes range within [1..M, 1..N].
'	M   -   number of rows in matrix A.
'	N   -   number of columns in matrix A.
'
'Output parameters:
'	A   -   matrices L and U in compact form (see below).
'			Array whose indexes range within [1..M, 1..N].
'	Pivots - permutation matrix in compact form (see below).
'			Array whose index ranges within [1..Min(M,N)].
'
'Matrix A is represented as A = P * L * U, where P is a permutation matrix,
'matrix L - lower triangular (or lower trapezoid, if M>N) matrix,
&#39;U - upper triangular &#40;or upper trapezoid, if M<N&#41; matrix.
&#39;
&#39;Let M be equal to 4 and N be equal to 3&#58;
&#39;
&#39;				   &#40;  1		  &#41;	&#40; U11 U12 U13  &#41;
&#39;A = P1 * P2 * P3 * &#40; L21  1	  &#41;  * &#40;	 U22 U23  &#41;
&#39;				   &#40; L31 L32  1  &#41;	&#40;		 U33  &#41;
&#39;				   &#40; L41 L42 L43 &#41;
&#39;
&#39;Matrix L has size MxMin&#40;M,N&#41;, matrix U has size Min&#40;M,N&#41;xN, matrix P&#40;i&#41; is
&#39;a permutation of the identity matrix of size MxM with numbers I and Pivots&#91;I&#93;.
&#39;
&#39;The algorithm returns array Pivots and the following matrix which replaces
&#39;matrix A and contains matrices L and U in compact form &#40;the example applies
&#39;to M=4, N=3&#41;.
&#39;
&#39; &#40; U11 U12 U13 &#41;
&#39; &#40; L21 U22 U23 &#41;
&#39; &#40; L31 L32 U33 &#41;
&#39; &#40; L41 L42 L43 &#41;
&#39;
&#39;As we can see, the unit diagonal isn&#39;t stored.
&#39;
&#39;  -- LAPACK routine &#40;version 3.0&#41; --
&#39;	 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
&#39;	 Courant Institute, Argonne National Lab, and Rice University
&#39;	 June 30, 1992
&#39;
&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
Public Sub LUDecomposition&#40;ByRef A&#40;&#41; As Double, _
		 ByVal M As Long, _
		 ByVal N As Long, _
		 ByRef Pivots&#40;&#41; As Long&#41;
	Dim I As Long
	Dim J As Long
	Dim JP As Long
	Dim T1&#40;&#41; As Double
	Dim s As Double
	Dim i_ As Long

	ReDim Pivots&#40;1# To MinInt&#40;M, N&#41;&#41;
	ReDim T1&#40;1# To MaxInt&#40;M, N&#41;&#41;
	
	&#39;
	&#39; Quick return if possible
	&#39;
	If M=0# Or N=0# then
		Exit Sub
	End If
	For J=1# To MinInt&#40;M, N&#41; Step 1
		
		&#39;
		&#39; Find pivot and test for singularity.
		&#39;
		JP = J
		For I=J+1# To M Step 1
			If Abs&#40;A&#40;I,J&#41;&#41;>Abs&#40;A&#40;JP,J&#41;&#41; then
				JP = I
			End If
		Next I
		Pivots&#40;J&#41; = JP
		If A&#40;JP,J&#41;<>0# then
			
			&#39;
			&#39;Apply the interchange to rows
			&#39;
			If JP<>J then
				For i_ = 1# To N Step 1
					T1&#40;i_&#41; = A&#40;J,i_&#41;
				Next i_
				For i_ = 1# To N Step 1
					A&#40;J,i_&#41; = A&#40;JP,i_&#41;
				Next i_
				For i_ = 1# To N Step 1
					A&#40;JP,i_&#41; = T1&#40;i_&#41;
				Next i_
			End If
			
			&#39;
			&#39;Compute elements J+1&#58;M of J-th column.
			&#39;
			If J<M then
				
				&#39;
				&#39; CALL DSCAL&#40; M-J, ONE / A&#40; J, J &#41;, A&#40; J+1, J &#41;, 1 &#41;
				&#39;
				JP = J+1#
				S = 1#/A&#40;J,J&#41;
				For i_ = JP To M Step 1
					A&#40;i_,J&#41; = S*A&#40;i_,J&#41;
				Next i_
			End If
		End If
		If J<MinInt&#40;M, N&#41; then
			
			&#39;
			&#39;Update trailing submatrix.
			&#39;CALL DGER&#40; M-J, N-J, -ONE, A&#40; J+1, J &#41;, 1, A&#40; J, J+1 &#41;, LDA,A&#40; J+1, J+1 &#41;, LDA &#41;
			&#39;
			JP = J+1#
			For I=J+1# To M Step 1
				S = A&#40;I,J&#41;
				For i_ = JP To N Step 1
					A&#40;I,i_&#41; = A&#40;I,i_&#41; - S*A&#40;J,i_&#41;
				Next i_
			Next I
		End If
	Next J
End Sub


&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
&#39;LU decomposition of a general matrix of size MxN
&#39;
&#39;It uses LUDecomposition. L and U are not output in compact form, but as
&#39;separate general matrices filled up by zero elements in their
&#39;corresponding positions.
&#39;
&#39;This subroutine described here only serves the purpose to show
&#39;how the result of ComplexLUDecomposition subroutine could be unpacked.
&#39;
&#39;  -- ALGLIB --
&#39;	 Copyright 2005 by Bochkanov Sergey
&#39;
&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
Public Sub LUDecompositionUnpacked&#40;ByRef A_&#40;&#41; As Double, _
		 ByVal M As Long, _
		 ByVal N As Long, _
		 ByRef L&#40;&#41; As Double, _
		 ByRef U&#40;&#41; As Double, _
		 ByRef Pivots&#40;&#41; As Long&#41;
	Dim A&#40;&#41; As Double
	Dim I As Long
	Dim J As Long
	Dim MinMN As Long
	A = A_

	If M=0# Or N=0# then
		Exit Sub
	End If
	MinMN = MinInt&#40;M, N&#41;
	ReDim L&#40;1# To M, 1# To MinMN&#41;
	ReDim U&#40;1# To MinMN, 1# To N&#41;
	Call LUDecomposition&#40;A, M, N, Pivots&#41;
	For I=1# To M Step 1
		For J=1# To MinMN Step 1
			If J>I then
				L&#40;I,J&#41; = 0#
			End If
			If J=I then
				L&#40;I,J&#41; = 1#
			End If
			If J<I then
				L&#40;I,J&#41; = A&#40;I,J&#41;
			End If
		Next J
	Next I
	For I=1# To MinMN Step 1
		For J=1# To N Step 1
			If J<I then
				U&#40;I,J&#41; = 0#
			End If
			If J>=I then
				U&#40;I,J&#41; = A&#40;I,J&#41;
			End If
		Next J
	Next I
End Sub

lu.cs

/***********************************************************************
This code is generated by the AlgoPascal translator

This code is distributed under the ALGLIB license
	&#40;see http&#58;//www.alglib.net/copyrules.php for details&#41;
***********************************************************************/
/*************************************************************************
LU decomposition of a general matrix of size MxN

The subroutine calculates the LU decomposition of a rectangular general
matrix with partial pivoting &#40;with row permutations&#41;.

Input parameters&#58;
	A   -   matrix A whose indexes range within &#91;1..M, 1..N&#93;.
	M   -   number of rows in matrix A.
	N   -   number of columns in matrix A.

Output parameters&#58;
	A   -   matrices L and U in compact form &#40;see below&#41;.
			Array whose indexes range within &#91;1..M, 1..N&#93;.
	Pivots - permutation matrix in compact form &#40;see below&#41;.
			Array whose index ranges within &#91;1..Min&#40;M,N&#41;&#93;.

Matrix A is represented as A = P * L * U, where P is a permutation matrix,
matrix L - lower triangular &#40;or lower trapezoid, if M>N&#41; matrix,
U - upper triangular &#40;or upper trapezoid, if M<N&#41; matrix.

Let M be equal to 4 and N be equal to 3&#58;

				   &#40;  1		  &#41;	&#40; U11 U12 U13  &#41;
A = P1 * P2 * P3 * &#40; L21  1	  &#41;  * &#40;	 U22 U23  &#41;
				   &#40; L31 L32  1  &#41;	&#40;		 U33  &#41;
				   &#40; L41 L42 L43 &#41;

Matrix L has size MxMin&#40;M,N&#41;, matrix U has size Min&#40;M,N&#41;xN, matrix P&#40;i&#41; is
a permutation of the identity matrix of size MxM with numbers I and Pivots&#91;I&#93;.

The algorithm returns array Pivots and the following matrix which replaces
matrix A and contains matrices L and U in compact form &#40;the example applies
to M=4, N=3&#41;.

 &#40; U11 U12 U13 &#41;
 &#40; L21 U22 U23 &#41;
 &#40; L31 L32 U33 &#41;
 &#40; L41 L42 L43 &#41;

As we can see, the unit diagonal isn&#39;t stored.

  -- LAPACK routine &#40;version 3.0&#41; --
	 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	 Courant Institute, Argonne National Lab, and Rice University
	 June 30, 1992
*************************************************************************/
public static void ludecomposition&#40;ref double&#91;,&#93; a,
	int m,
	int n,
	ref int&#91;&#93; pivots&#41;
{
	int i = 0;
	int j = 0;
	int jp = 0;
	double&#91;&#93; t1 = new double&#91;0&#93;;
	double s = 0;
	int i_ = 0;


	pivots = new int&#91;Math.Min&#40;m, n&#41;+1&#93;;
	t1 = new double&#91;Math.Max&#40;m, n&#41;+1&#93;;
	System.Diagnostics.Debug.Assert&#40;m>=0 & n>=0, &#34;Error in LUDecomposition&#58; incorrect function arguments&#34;&#41;;
	
	//
	// Quick return if possible
	//
	if&#40; m==0 | n==0 &#41;
	{
		return;
	}
	for&#40;j=1; j<=Math.Min&#40;m, n&#41;; j++&#41;
	{
		
		//
		// Find pivot and test for singularity.
		//
		jp = j;
		for&#40;i=j+1; i<=m; i++&#41;
		{
			if&#40; Math.Abs&#40;a&#91;i,j&#93;&#41;>Math.Abs&#40;a&#91;jp,j&#93;&#41; &#41;
			{
				jp = i;
			}
		}
		pivots&#91;j&#93; = jp;
		if&#40; a&#91;jp,j&#93;!=0 &#41;
		{
			
			//
			//Apply the interchange to rows
			//
			if&#40; jp!=j &#41;
			{
				for&#40;i_=1; i_<=n;i_++&#41;
				{
					t1&#91;i_&#93; = a&#91;j,i_&#93;;
				}
				for&#40;i_=1; i_<=n;i_++&#41;
				{
					a&#91;j,i_&#93; = a&#91;jp,i_&#93;;
				}
				for&#40;i_=1; i_<=n;i_++&#41;
				{
					a&#91;jp,i_&#93; = t1&#91;i_&#93;;
				}
			}
			
			//
			//Compute elements J+1&#58;M of J-th column.
			//
			if&#40; j<m &#41;
			{
				
				//
				// CALL DSCAL&#40; M-J, ONE / A&#40; J, J &#41;, A&#40; J+1, J &#41;, 1 &#41;
				//
				jp = j+1;
				s = 1/a&#91;j,j&#93;;
				for&#40;i_=jp; i_<=m;i_++&#41;
				{
					a&#91;i_,j&#93; = s*a&#91;i_,j&#93;;
				}
			}
		}
		if&#40; j<Math.Min&#40;m, n&#41; &#41;
		{
			
			//
			//Update trailing submatrix.
			//CALL DGER&#40; M-J, N-J, -ONE, A&#40; J+1, J &#41;, 1, A&#40; J, J+1 &#41;, LDA,A&#40; J+1, J+1 &#41;, LDA &#41;
			//
			jp = j+1;
			for&#40;i=j+1; i<=m; i++&#41;
			{
				s = a&#91;i,j&#93;;
				for&#40;i_=jp; i_<=n;i_++&#41;
				{
					a&#91;i,i_&#93; = a&#91;i,i_&#93; - s*a&#91;j,i_&#93;;
				}
			}
		}
	}
}


/*************************************************************************
LU decomposition of a general matrix of size MxN

It uses LUDecomposition. L and U are not output in compact form, but as
separate general matrices filled up by zero elements in their
corresponding positions.

This subroutine described here only serves the purpose to show
how the result of ComplexLUDecomposition subroutine could be unpacked.

  -- ALGLIB --
	 Copyright 2005 by Bochkanov Sergey
*************************************************************************/
public static void ludecompositionunpacked&#40;double&#91;,&#93; a,
	int m,
	int n,
	ref double&#91;,&#93; l,
	ref double&#91;,&#93; u,
	ref int&#91;&#93; pivots&#41;
{
	int i = 0;
	int j = 0;
	int minmn = 0;


	a = &#40;double&#91;,&#93;&#41;a.Clone&#40;&#41;;


	if&#40; m==0 | n==0 &#41;
	{
		return;
	}
	minmn = Math.Min&#40;m, n&#41;;
	l = new double&#91;m+1, minmn+1&#93;;
	u = new double&#91;minmn+1, n+1&#93;;
	ludecomposition&#40;ref a, m, n, ref pivots&#41;;
	for&#40;i=1; i<=m; i++&#41;
	{
		for&#40;j=1; j<=minmn; j++&#41;
		{
			if&#40; j>i &#41;
			{
				l&#91;i,j&#93; = 0;
			}
			if&#40; j==i &#41;
			{
				l&#91;i,j&#93; = 1;
			}
			if&#40; j<i &#41;
			{
				l&#91;i,j&#93; = a&#91;i,j&#93;;
			}
		}
	}
	for&#40;i=1; i<=minmn; i++&#41;
	{
		for&#40;j=1; j<=n; j++&#41;
		{
			if&#40; j<i &#41;
			{
				u&#91;i,j&#93; = 0;
			}
			if&#40; j>=i &#41;
			{
				u&#91;i,j&#93; = a&#91;i,j&#93;;
			}
		}
	}
}

lu.cpp

/***********************************************************************
This code is generated by the AlgoPascal translator

This code is distributed under the ALGLIB license
	&#40;see http&#58;//www.alglib.net/copyrules.php for details&#41;
***********************************************************************/

#include &#34;ap.h&#34;

void ludecomposition&#40;ap&#58;&#58;real_2d_array& a,
	 int m,
	 int n,
	 ap&#58;&#58;integer_1d_array& pivots&#41;;
void ludecompositionunpacked&#40;ap&#58;&#58;real_2d_array a,
	 int m,
	 int n,
	 ap&#58;&#58;real_2d_array& l,
	 ap&#58;&#58;real_2d_array& u,
	 ap&#58;&#58;integer_1d_array& pivots&#41;;

/*************************************************************************
LU decomposition of a general matrix of size MxN

The subroutine calculates the LU decomposition of a rectangular general
matrix with partial pivoting &#40;with row permutations&#41;.

Input parameters&#58;
	A   -   matrix A whose indexes range within &#91;1..M, 1..N&#93;.
	M   -   number of rows in matrix A.
	N   -   number of columns in matrix A.

Output parameters&#58;
	A   -   matrices L and U in compact form &#40;see below&#41;.
			Array whose indexes range within &#91;1..M, 1..N&#93;.
	Pivots - permutation matrix in compact form &#40;see below&#41;.
			Array whose index ranges within &#91;1..Min&#40;M,N&#41;&#93;.

Matrix A is represented as A = P * L * U, where P is a permutation matrix,
matrix L - lower triangular &#40;or lower trapezoid, if M>N&#41; matrix,
U - upper triangular &#40;or upper trapezoid, if M<N&#41; matrix.

Let M be equal to 4 and N be equal to 3&#58;

				   &#40;  1		  &#41;	&#40; U11 U12 U13  &#41;
A = P1 * P2 * P3 * &#40; L21  1	  &#41;  * &#40;	 U22 U23  &#41;
				   &#40; L31 L32  1  &#41;	&#40;		 U33  &#41;
				   &#40; L41 L42 L43 &#41;

Matrix L has size MxMin&#40;M,N&#41;, matrix U has size Min&#40;M,N&#41;xN, matrix P&#40;i&#41; is
a permutation of the identity matrix of size MxM with numbers I and Pivots&#91;I&#93;.

The algorithm returns array Pivots and the following matrix which replaces
matrix A and contains matrices L and U in compact form &#40;the example applies
to M=4, N=3&#41;.

 &#40; U11 U12 U13 &#41;
 &#40; L21 U22 U23 &#41;
 &#40; L31 L32 U33 &#41;
 &#40; L41 L42 L43 &#41;

As we can see, the unit diagonal isn&#39;t stored.

  -- LAPACK routine &#40;version 3.0&#41; --
	 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	 Courant Institute, Argonne National Lab, and Rice University
	 June 30, 1992
*************************************************************************/
void ludecomposition&#40;ap&#58;&#58;real_2d_array& a,
	 int m,
	 int n,
	 ap&#58;&#58;integer_1d_array& pivots&#41;
{
	int i;
	int j;
	int jp;
	ap&#58;&#58;real_1d_array t1;
	double s;

	pivots.setbounds&#40;1, ap&#58;&#58;minint&#40;m, n&#41;&#41;;
	t1.setbounds&#40;1, ap&#58;&#58;maxint&#40;m, n&#41;&#41;;
	ap&#58;&#58;ap_error&#58;&#58;make_assertion&#40;m>=0&&n>=0&#41;;
	
	//
	// Quick return if possible
	//
	if&#40; m==0||n==0 &#41;
	{
		return;
	}
	for&#40;j = 1; j <= ap&#58;&#58;minint&#40;m, n&#41;; j++&#41;
	{
		
		//
		// Find pivot and test for singularity.
		//
		jp = j;
		for&#40;i = j+1; i <= m; i++&#41;
		{
			if&#40; fabs&#40;a&#40;i,j&#41;&#41;>fabs&#40;a&#40;jp,j&#41;&#41; &#41;
			{
				jp = i;
			}
		}
		pivots&#40;j&#41; = jp;
		if&#40; a&#40;jp,j&#41;!=0 &#41;
		{
			
			//
			//Apply the interchange to rows
			//
			if&#40; jp!=j &#41;
			{
				ap&#58;&#58;vmove&#40;t1.getvector&#40;1, n&#41;, a.getrow&#40;j, 1, n&#41;&#41;;
				ap&#58;&#58;vmove&#40;a.getrow&#40;j, 1, n&#41;, a.getrow&#40;jp, 1, n&#41;&#41;;
				ap&#58;&#58;vmove&#40;a.getrow&#40;jp, 1, n&#41;, t1.getvector&#40;1, n&#41;&#41;;
			}
			
			//
			//Compute elements J+1&#58;M of J-th column.
			//
			if&#40; j<m &#41;
			{
				
				//
				// CALL DSCAL&#40; M-J, ONE / A&#40; J, J &#41;, A&#40; J+1, J &#41;, 1 &#41;
				//
				jp = j+1;
				s = 1/a&#40;j,j&#41;;
				ap&#58;&#58;vmul&#40;a.getcolumn&#40;j, jp, m&#41;, s&#41;;
			}
		}
		if&#40; j<ap&#58;&#58;minint&#40;m, n&#41; &#41;
		{
			
			//
			//Update trailing submatrix.
			//CALL DGER&#40; M-J, N-J, -ONE, A&#40; J+1, J &#41;, 1, A&#40; J, J+1 &#41;, LDA,A&#40; J+1, J+1 &#41;, LDA &#41;
			//
			jp = j+1;
			for&#40;i = j+1; i <= m; i++&#41;
			{
				s = a&#40;i,j&#41;;
				ap&#58;&#58;vsub&#40;a.getrow&#40;i, jp, n&#41;, a.getrow&#40;j, jp, n&#41;, s&#41;;
			}
		}
	}
}


/*************************************************************************
LU decomposition of a general matrix of size MxN

It uses LUDecomposition. L and U are not output in compact form, but as
separate general matrices filled up by zero elements in their
corresponding positions.

This subroutine described here only serves the purpose to show
how the result of ComplexLUDecomposition subroutine could be unpacked.

  -- ALGLIB --
	 Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void ludecompositionunpacked&#40;ap&#58;&#58;real_2d_array a,
	 int m,
	 int n,
	 ap&#58;&#58;real_2d_array& l,
	 ap&#58;&#58;real_2d_array& u,
	 ap&#58;&#58;integer_1d_array& pivots&#41;
{
	int i;
	int j;
	int minmn;

	if&#40; m==0||n==0 &#41;
	{
		return;
	}
	minmn = ap&#58;&#58;minint&#40;m, n&#41;;
	l.setbounds&#40;1, m, 1, minmn&#41;;
	u.setbounds&#40;1, minmn, 1, n&#41;;
	ludecomposition&#40;a, m, n, pivots&#41;;
	for&#40;i = 1; i <= m; i++&#41;
	{
		for&#40;j = 1; j <= minmn; j++&#41;
		{
			if&#40; j>i &#41;
			{
				l&#40;i,j&#41; = 0;
			}
			if&#40; j==i &#41;
			{
				l&#40;i,j&#41; = 1;
			}
			if&#40; j<i &#41;
			{
				l&#40;i,j&#41; = a&#40;i,j&#41;;
			}
		}
	}
	for&#40;i = 1; i <= minmn; i++&#41;
	{
		for&#40;j = 1; j <= n; j++&#41;
		{
			if&#40; j<i &#41;
			{
				u&#40;i,j&#41; = 0;
			}
			if&#40; j>=i &#41;
			{
				u&#40;i,j&#41; = a&#40;i,j&#41;;
			}
		}
	}
}

Anh em cố mà nghiên cứu nhá :">

Trong chủ đề: Nghịch đảo ma trận bằng các ngôn ngữ

24-08-2007 - 00:57

invtriangular.bas

&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
&#39;This code is generated by the AlgoPascal translator
&#39;
&#39;This code is distributed under the ALGLIB license
&#39;	&#40;see http&#58;//www.alglib.net/copyrules.php for details&#41;
&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
&#39;Routines
&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
&#39;Triangular matrix inversion
&#39;
&#39;The subroutine inverts the following types of matrices&#58;
&#39;	* upper triangular
&#39;	* upper triangular with unit diagonal
&#39;	* lower triangular
&#39;	* lower triangular with unit diagonal
&#39;
&#39;In case of an upper &#40;lower&#41; triangular matrix,  the  inverse  matrix  will
&#39;also be upper &#40;lower&#41; triangular, and after the end of the algorithm,  the
&#39;inverse matrix replaces the source matrix. The elements  below &#40;above&#41; the
&#39;main diagonal are not changed by the algorithm.
&#39;
&#39;If  the matrix  has a unit diagonal, the inverse matrix also  has  a  unit
&#39;diagonal, and the diagonal elements are not passed to the algorithm.
&#39;
&#39;Input parameters&#58;
&#39;	A	   -   matrix.
&#39;				Array whose indexes range within &#91;1..N, 1..N&#93;.
&#39;	N	   -   size of matrix A.
&#39;	IsUpper -   True, if the matrix is upper triangular.
&#39;	IsUnitTriangular
&#39;			-   True, if the matrix has a unit diagonal.
&#39;
&#39;Output parameters&#58;
&#39;	A	   -   inverse matrix &#40;if the problem is not degenerate&#41;.
&#39;
&#39;Result&#58;
&#39;	True, if the matrix is not singular.
&#39;	False, if the matrix is singular.
&#39;
&#39;  -- LAPACK routine &#40;version 3.0&#41; --
&#39;	 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
&#39;	 Courant Institute, Argonne National Lab, and Rice University
&#39;	 February 29, 1992
&#39;
&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;&#39;
Public Function InvTriangular&#40;ByRef A&#40;&#41; As Double, _
		 ByVal N As Long, _
		 ByVal IsUpper As Boolean, _
		 ByVal IsUnitTriangular As Boolean&#41; As Boolean
	Dim Result As Boolean
	Dim NOUNIT As Boolean
	Dim I As Long
	Dim J As Long
	Dim NMJ As Long
	Dim JM1 As Long
	Dim JP1 As Long
	Dim V As Double
	Dim AJJ As Double
	Dim T&#40;&#41; As Double
	Dim i_ As Long

	Result = True
	ReDim T&#40;1# To N&#41;
	
	&#39;
	&#39; Test the input parameters.
	&#39;
	NOUNIT =  Not IsUnitTriangular
	If IsUpper then
		
		&#39;
		&#39; Compute inverse of upper triangular matrix.
		&#39;
		For J=1# To N Step 1
			If NOUNIT then
				If A&#40;J,J&#41;=0# then
					Result = False
					InvTriangular = Result
					Exit Function
				End If
				A&#40;J,J&#41; = 1#/A&#40;J,J&#41;
				AJJ = -A&#40;J,J&#41;
			Else
				AJJ = -1#
			End If
			
			&#39;
			&#39; Compute elements 1&#58;j-1 of j-th column.
			&#39;
			If J>1# then
				JM1 = J-1#
				For i_ = 1# To JM1 Step 1
					T&#40;i_&#41; = A&#40;i_,J&#41;
				Next i_
				For I=1# To J-1# Step 1
					If I<J-1# then
						V = 0.0
						For i_ = I+1# To JM1 Step 1
							V = V + A&#40;I,i_&#41;*T&#40;i_&#41;
						Next i_
					Else
						V = 0#
					End If
					If NOUNIT then
						A&#40;I,J&#41; = V+A&#40;I,I&#41;*T&#40;I&#41;
					Else
						A&#40;I,J&#41; = V+T&#40;I&#41;
					End If
				Next I
				For i_ = 1# To JM1 Step 1
					A&#40;i_,J&#41; = AJJ*A&#40;i_,J&#41;
				Next i_
			End If
		Next J
	Else
		
		&#39;
		&#39; Compute inverse of lower triangular matrix.
		&#39;
		For J=N To 1# Step -1
			If NOUNIT then
				If A&#40;J,J&#41;=0# then
					Result = False
					InvTriangular = Result
					Exit Function
				End If
				A&#40;J,J&#41; = 1#/A&#40;J,J&#41;
				AJJ = -A&#40;J,J&#41;
			Else
				AJJ = -1#
			End If
			If J<N then
				
				&#39;
				&#39; Compute elements j+1&#58;n of j-th column.
				&#39;
				NMJ = N-J
				JP1 = J+1#
				For i_ = JP1 To N Step 1
					T&#40;i_&#41; = A&#40;i_,J&#41;
				Next i_
				For I=J+1# To N Step 1
					If I>J+1# then
						V = 0.0
						For i_ = JP1 To I-1# Step 1
							V = V + A&#40;I,i_&#41;*T&#40;i_&#41;
						Next i_
					Else
						V = 0#
					End If
					If NOUNIT then
						A&#40;I,J&#41; = V+A&#40;I,I&#41;*T&#40;I&#41;
					Else
						A&#40;I,J&#41; = V+T&#40;I&#41;
					End If
				Next I
				For i_ = JP1 To N Step 1
					A&#40;i_,J&#41; = AJJ*A&#40;i_,J&#41;
				Next i_
			End If
		Next J
	End If

	InvTriangular = Result
End Function

invtriangular.cs

/***********************************************************************
This code is generated by the AlgoPascal translator

This code is distributed under the ALGLIB license
	&#40;see http&#58;//www.alglib.net/copyrules.php for details&#41;
***********************************************************************/
/*************************************************************************
Triangular matrix inversion

The subroutine inverts the following types of matrices&#58;
	* upper triangular
	* upper triangular with unit diagonal
	* lower triangular
	* lower triangular with unit diagonal

In case of an upper &#40;lower&#41; triangular matrix,  the  inverse  matrix  will
also be upper &#40;lower&#41; triangular, and after the end of the algorithm,  the
inverse matrix replaces the source matrix. The elements  below &#40;above&#41; the
main diagonal are not changed by the algorithm.

If  the matrix  has a unit diagonal, the inverse matrix also  has  a  unit
diagonal, and the diagonal elements are not passed to the algorithm.

Input parameters&#58;
	A	   -   matrix.
				Array whose indexes range within &#91;1..N, 1..N&#93;.
	N	   -   size of matrix A.
	IsUpper -   True, if the matrix is upper triangular.
	IsUnitTriangular
			-   True, if the matrix has a unit diagonal.

Output parameters&#58;
	A	   -   inverse matrix &#40;if the problem is not degenerate&#41;.

Result&#58;
	True, if the matrix is not singular.
	False, if the matrix is singular.

  -- LAPACK routine &#40;version 3.0&#41; --
	 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	 Courant Institute, Argonne National Lab, and Rice University
	 February 29, 1992
*************************************************************************/
public static bool invtriangular&#40;ref double&#91;,&#93; a,
	int n,
	bool isupper,
	bool isunittriangular&#41;
{
	bool result = new bool&#40;&#41;;
	bool nounit = new bool&#40;&#41;;
	int i = 0;
	int j = 0;
	int nmj = 0;
	int jm1 = 0;
	int jp1 = 0;
	double v = 0;
	double ajj = 0;
	double&#91;&#93; t = new double&#91;0&#93;;
	int i_ = 0;


	result = true;
	t = new double&#91;n+1&#93;;
	
	//
	// Test the input parameters.
	//
	nounit = !isunittriangular;
	if&#40; isupper &#41;
	{
		
		//
		// Compute inverse of upper triangular matrix.
		//
		for&#40;j=1; j<=n; j++&#41;
		{
			if&#40; nounit &#41;
			{
				if&#40; a&#91;j,j&#93;==0 &#41;
				{
					result = false;
					return result;
				}
				a&#91;j,j&#93; = 1/a&#91;j,j&#93;;
				ajj = -a&#91;j,j&#93;;
			}
			else
			{
				ajj = -1;
			}
			
			//
			// Compute elements 1&#58;j-1 of j-th column.
			//
			if&#40; j>1 &#41;
			{
				jm1 = j-1;
				for&#40;i_=1; i_<=jm1;i_++&#41;
				{
					t&#91;i_&#93; = a&#91;i_,j&#93;;
				}
				for&#40;i=1; i<=j-1; i++&#41;
				{
					if&#40; i<j-1 &#41;
					{
						v = 0.0;
						for&#40;i_=i+1; i_<=jm1;i_++&#41;
						{
							v += a&#91;i,i_&#93;*t&#91;i_&#93;;
						}
					}
					else
					{
						v = 0;
					}
					if&#40; nounit &#41;
					{
						a&#91;i,j&#93; = v+a&#91;i,i&#93;*t&#91;i&#93;;
					}
					else
					{
						a&#91;i,j&#93; = v+t&#91;i&#93;;
					}
				}
				for&#40;i_=1; i_<=jm1;i_++&#41;
				{
					a&#91;i_,j&#93; = ajj*a&#91;i_,j&#93;;
				}
			}
		}
	}
	else
	{
		
		//
		// Compute inverse of lower triangular matrix.
		//
		for&#40;j=n; j>=1; j--&#41;
		{
			if&#40; nounit &#41;
			{
				if&#40; a&#91;j,j&#93;==0 &#41;
				{
					result = false;
					return result;
				}
				a&#91;j,j&#93; = 1/a&#91;j,j&#93;;
				ajj = -a&#91;j,j&#93;;
			}
			else
			{
				ajj = -1;
			}
			if&#40; j<n &#41;
			{
				
				//
				// Compute elements j+1&#58;n of j-th column.
				//
				nmj = n-j;
				jp1 = j+1;
				for&#40;i_=jp1; i_<=n;i_++&#41;
				{
					t&#91;i_&#93; = a&#91;i_,j&#93;;
				}
				for&#40;i=j+1; i<=n; i++&#41;
				{
					if&#40; i>j+1 &#41;
					{
						v = 0.0;
						for&#40;i_=jp1; i_<=i-1;i_++&#41;
						{
							v += a&#91;i,i_&#93;*t&#91;i_&#93;;
						}
					}
					else
					{
						v = 0;
					}
					if&#40; nounit &#41;
					{
						a&#91;i,j&#93; = v+a&#91;i,i&#93;*t&#91;i&#93;;
					}
					else
					{
						a&#91;i,j&#93; = v+t&#91;i&#93;;
					}
				}
				for&#40;i_=jp1; i_<=n;i_++&#41;
				{
					a&#91;i_,j&#93; = ajj*a&#91;i_,j&#93;;
				}
			}
		}
	}
	return result;
}

invtriangular.cpp

/***********************************************************************
This code is generated by the AlgoPascal translator

This code is distributed under the ALGLIB license
	&#40;see http&#58;//www.alglib.net/copyrules.php for details&#41;
***********************************************************************/

#include &#34;ap.h&#34;

bool invtriangular&#40;ap&#58;&#58;real_2d_array& a,
	 int n,
	 bool isupper,
	 bool isunittriangular&#41;;

/*************************************************************************
Triangular matrix inversion

The subroutine inverts the following types of matrices&#58;
	* upper triangular
	* upper triangular with unit diagonal
	* lower triangular
	* lower triangular with unit diagonal

In case of an upper &#40;lower&#41; triangular matrix,  the  inverse  matrix  will
also be upper &#40;lower&#41; triangular, and after the end of the algorithm,  the
inverse matrix replaces the source matrix. The elements  below &#40;above&#41; the
main diagonal are not changed by the algorithm.

If  the matrix  has a unit diagonal, the inverse matrix also  has  a  unit
diagonal, and the diagonal elements are not passed to the algorithm.

Input parameters&#58;
	A	   -   matrix.
				Array whose indexes range within &#91;1..N, 1..N&#93;.
	N	   -   size of matrix A.
	IsUpper -   True, if the matrix is upper triangular.
	IsUnitTriangular
			-   True, if the matrix has a unit diagonal.

Output parameters&#58;
	A	   -   inverse matrix &#40;if the problem is not degenerate&#41;.

Result&#58;
	True, if the matrix is not singular.
	False, if the matrix is singular.

  -- LAPACK routine &#40;version 3.0&#41; --
	 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	 Courant Institute, Argonne National Lab, and Rice University
	 February 29, 1992
*************************************************************************/
bool invtriangular&#40;ap&#58;&#58;real_2d_array& a,
	 int n,
	 bool isupper,
	 bool isunittriangular&#41;
{
	bool result;
	bool nounit;
	int i;
	int j;
	int nmj;
	int jm1;
	int jp1;
	double v;
	double ajj;
	ap&#58;&#58;real_1d_array t;

	result = true;
	t.setbounds&#40;1, n&#41;;
	
	//
	// Test the input parameters.
	//
	nounit = !isunittriangular;
	if&#40; isupper &#41;
	{
		
		//
		// Compute inverse of upper triangular matrix.
		//
		for&#40;j = 1; j <= n; j++&#41;
		{
			if&#40; nounit &#41;
			{
				if&#40; a&#40;j,j&#41;==0 &#41;
				{
					result = false;
					return result;
				}
				a&#40;j,j&#41; = 1/a&#40;j,j&#41;;
				ajj = -a&#40;j,j&#41;;
			}
			else
			{
				ajj = -1;
			}
			
			//
			// Compute elements 1&#58;j-1 of j-th column.
			//
			if&#40; j>1 &#41;
			{
				jm1 = j-1;
				ap&#58;&#58;vmove&#40;t.getvector&#40;1, jm1&#41;, a.getcolumn&#40;j, 1, jm1&#41;&#41;;
				for&#40;i = 1; i <= j-1; i++&#41;
				{
					if&#40; i<j-1 &#41;
					{
						v = ap&#58;&#58;vdotproduct&#40;a.getrow&#40;i, i+1, jm1&#41;, t.getvector&#40;i+1, jm1&#41;&#41;;
					}
					else
					{
						v = 0;
					}
					if&#40; nounit &#41;
					{
						a&#40;i,j&#41; = v+a&#40;i,i&#41;*t&#40;i&#41;;
					}
					else
					{
						a&#40;i,j&#41; = v+t&#40;i&#41;;
					}
				}
				ap&#58;&#58;vmul&#40;a.getcolumn&#40;j, 1, jm1&#41;, ajj&#41;;
			}
		}
	}
	else
	{
		
		//
		// Compute inverse of lower triangular matrix.
		//
		for&#40;j = n; j >= 1; j--&#41;
		{
			if&#40; nounit &#41;
			{
				if&#40; a&#40;j,j&#41;==0 &#41;
				{
					result = false;
					return result;
				}
				a&#40;j,j&#41; = 1/a&#40;j,j&#41;;
				ajj = -a&#40;j,j&#41;;
			}
			else
			{
				ajj = -1;
			}
			if&#40; j<n &#41;
			{
				
				//
				// Compute elements j+1&#58;n of j-th column.
				//
				nmj = n-j;
				jp1 = j+1;
				ap&#58;&#58;vmove&#40;t.getvector&#40;jp1, n&#41;, a.getcolumn&#40;j, jp1, n&#41;&#41;;
				for&#40;i = j+1; i <= n; i++&#41;
				{
					if&#40; i>j+1 &#41;
					{
						v = ap&#58;&#58;vdotproduct&#40;a.getrow&#40;i, jp1, i-1&#41;, t.getvector&#40;jp1, i-1&#41;&#41;;
					}
					else
					{
						v = 0;
					}
					if&#40; nounit &#41;
					{
						a&#40;i,j&#41; = v+a&#40;i,i&#41;*t&#40;i&#41;;
					}
					else
					{
						a&#40;i,j&#41; = v+t&#40;i&#41;;
					}
				}
				ap&#58;&#58;vmul&#40;a.getcolumn&#40;j, jp1, n&#41;, ajj&#41;;
			}
		}
	}
	return result;
}

  1. Diễn đàn Toán học
  2. → Đang xem trang cá nhân: Bài viết: sans_amour