محاسبه دترمینان به روش گاوس
اینم مثل قبلی
کد:
/****************************************************
* Calculate the determinant of a real square matrix *
* A(n,n) by Gauss method with full pivoting. *
* ------------------------------------------------- *
* Ref.: "Algèbre - Algorithmes et programmes en *
* Pascal By Jean-Louis Jardrin, Dunod - *
* Bordas Paris, 1988 p. 76-79" [BIBLI 10]. *
* ------------------------------------------------- *
* SAMPLE RUN: *
* (Calculate the determinant of matrix: *
* 10 18 1 14 22 *
* 4 12 25 8 16 *
* 23 6 19 2 15 *
* 17 5 13 21 9 *
* 11 24 7 20 3 ) *
* *
* Input size of square real matrix: 5 *
* *
* Line 1 *
* Element 1: 10 *
* Element 2: 18 *
* Element 3: 1 *
* Element 4: 14 *
* Element 5: 22 *
* *
* Line 2 *
* Element 1: 4 *
* Element 2: 12 *
* Element 3: 25 *
* Element 4: 8 *
* Element 5: 16 *
* *
* Line 3 *
* Element 1: 23 *
* Element 2: 6 *
* Element 3: 19 *
* Element 4: 2 *
* Element 5: 15 *
* *
* Line 4 *
* Element 1: 17 *
* Element 2: 5 *
* Element 3: 13 *
* Element 4: 21 *
* Element 5: 9 *
* *
* Line 5 *
* Element 1: 11 *
* Element 2: 24 *
* Element 3: 7 *
* Element 4: 20 *
* Element 5: 3 *
* *
* Determinant = -4680000.000000 *
* *
* C++ version with dynamic allocations *
* by Jean-Pierre Moreau, Paris. *
* *
* (To link with basis_r.cpp and vmblock.cpp). *
****************************************************/
#include <basis.h>
#include <vmblock.h>
int n; //size of matrix A
REAL **A; //pointer to input matrix
REAL det; //determinant of matrix A
REAL eps; //desired precision
void *vmblock; //for dynamic allocation of A(n,n)
/*The function TSRGT applies to input real square matrix A(n,n) the upper
triangularization algorithm of Gauss method with full pivoting and keeps
trace of successive transformations done in integer vectors KP and LP.
----------------------------------------------------------------------------
Input parameters:
eps precision (real)
n size of A matrix (integer)
A pointer to input real square matrix (REAL**)
Output parameters:
it flag=1 if A matrix ok, =0 if A matrix is singular (integer)
C pointer to table storing main diagonal elements and supra-
diagonal elements of upper triangular matrix and the multi-
plying coefficients used during triangularization process (REAL**)
KP table storing informations concerning the column exchanges
during process (int*)
LP table storing informations concerning the line exchanges
during process (int*)
-----------------------------------------------------------------------------
The table C is first initialized to A matrix, then receives at each step k
of the triangularization process, usefull elements of A matrix at step k for
k=1,2,...n.
The variables po(real), lo and ko(integer) store respectively pivot at step k,
its line number and its column number.
-----------------------------------------------------------------------------*/
void TSRGT(REAL eps,int n,REAL **A,int *it,REAL **C, int *Kp,int *Lp) **
int i,j,k,ko,lo; REAL po,t0;
//set matrix C to A
for (i=1; i<=n; i++)
for (j=1; j<=n; j++)
C[i][j] = A[i][j];
*it=1; k=1;
while (*it==1 && k<n) **
po=C[k][k]; lo=k; ko=k;
for (i=k; i<=n; i++)
for (j=k; j<=n; j++)
if (fabs(C[i][j])>fabs(po)) **
po=C[i][j]; lo=i; ko=j;
}
Lp[k]=lo; Kp[k]=ko;
if (fabs(po)<eps) **
*it=0;
printf("\n Error: Pivot too small !\n");
}
else **
if (lo!=k)
for (j=k; j<=n; j++) **
t0=C[k][j]; C[k][j]=C[lo][j]; C[lo][j]=t0;
}
if (ko!=k)
for (i=1; i<=n; i++) **
t0=C[i][k]; C[i][k]=C[i][ko]; C[i][ko]=t0;
}
for (i=k+1; i<=n; i++) **
C[i][k] /= po;
for (j=k+1; j<=n; j++)
C[i][j] -= C[i][k]*C[k][j];
}
k++;
}
}
if (*it==1 && fabs(C[n][n])<eps) *it=0;
} //TSRGT()
/*The function DMGT returns the determinant of a real square matrix
A(n,n) by Gauss method with full pivoting.
----------------------------------------------------------------------------
Input parameters:
eps precision (REAL)
n size of A matrix (int)
A pointer to input real square matrix (REAL**)
Output parameters:
None
-----------------------------------------------------------------------------
The procedure TSRGT is used to reduce A matrix to an upper triangular matrix.
Output variables are it(integer), C(pMat), Kp and Lp(pVecI).
If it=0, matrix A is singular, if it=1, matrix A is regular. Table C contains
at location i,j (j>=i) the corresponding element of the upper triangular matrix.
Tables Lp and Kp contain informations relative to exchanges of line or column
that occured during the process. For instance, the element number k of Lp is
an integer <> k if an exchange of line has been made at step k (k=1,2,...,n).
The number of exchanges of lines and columns is stored in L(integer). the
determinant of A matrix is stored in d0(real).
Note: types pMat and pVecI must be declared and allocated for by main program,
except local variables C,Kp,Lp allocated (and disposed of) here.
-----------------------------------------------------------------------------*/
REAL DMGT(REAL eps,int n,REAL **A) **
int it,k,l;
REAL d0;
REAL **C;
int *Kp, *Lp;
void *vmblock1;
//allocate local matrix and vectors
vmblock1 = vminit();
C = (REAL **)vmalloc(vmblock1, MATRIX, n+1, n+1);
Kp = (int *) vmalloc(vmblock1, VEKTOR, n+1, 0);
Lp = (int *) vmalloc(vmblock1, VEKTOR, n+1, 0);
if (! vmcomplete(vmblock1))
**
LogError ("Memory full !", 0, __FILE__, __LINE__);
return 0;
}
TSRGT(eps,n,A,&it,C,Kp,Lp); //call triangularization procedure
if (it==0)
return 0.0; //matrix singular, det=0
else ** //matrix regular, det<>0
d0=1.0;
for (k=1; k<=n; k++) d0 *= C[k][k];
l=0;
for (k=1; k<n; k++) **
if (Lp[k]!=k) l++;
if (Kp[k]!=k) l++;
}
if ((l%2)!=0) d0=-d0;
}
return d0; //return determinant
vmfree(vmblock1);
}
void main() **
int i,j; REAL temp;
printf("\n Calculate the determinant of a real square matrix");
printf("\n by Gauss method with full pivoting.\n");
printf("\n Input size of square real matrix: "); scanf("%d", &n);
printf("\n");
//fix precision
eps=1e-10;
//allocate memory space for input matrix
vmblock = vminit();
A = (REAL **) vmalloc(vmblock, MATRIX, n+1, n+1);
if (! vmcomplete(vmblock))
**
LogError ("Memory full !", 0, __FILE__, __LINE__);
return;
}
for (i=1; i<=n; i++) **
printf(" Line %d\n",i);
for (j=1; j<=n; j++) **
printf(" Element %d: ", j); scanf("%lf", &temp);
A[i][j]=temp;
}
printf("\n");
}
det=DMGT(eps,n,A);
printf("\n Determinant = %f\n\n\n", det);
vmfree(vmblock);
}
//End of file deter.cpp
برنامه Cpp برای محاسبه دترمینان n*n به روش تجزیه LU
این رو هم من امتحان نکرده ام.
کد:
/****************************************************
* Calculate the determinant of a real square matrix *
* by the LU decomposition method *
* ------------------------------------------------- *
* SAMPLE RUN: *
* (Calculate the determinant of real square matrix: *
* | 1 0.5 0 0 0 0 0 0 0 | *
* |-1 0.5 1 0 0 0 0 0 0 | *
* | 0 -1 0 2.5 0 0 0 0 0 | *
* | 0 0 -1 -1.5 4 0 0 0 0 | *
* | 0 0 0 -1 -3 5.5 0 0 0 | *
* | 0 0 0 0 -1 -4.5 7 0 0 | *
* | 0 0 0 0 0 -1 -6 8.5 0 | *
* | 0 0 0 0 0 0 -1 -7.5 10 | *
* | 0 0 0 0 0 0 0 -1 -9 | ) *
* *
* *
* Determinant = 1.000000 *
* *
* *
* C++ version by J-P Moreau, Paris. *
*****************************************************
Exact value is: 1
---------------------------------------------------*/
#include <basis.h>
#include <vmblock.h>
#include "lu.h"
int main() **
REAL **A; // input matrix of size n x n
int *INDX; // dummy integer vector
void *vmblock = NULL;
// NOTA: index zero not used here.
int i, id, j, n, rc;
REAL det;
n=9;
vmblock = vminit();
A = (REAL **)vmalloc(vmblock, MATRIX, n+1, n+1);
INDX = (int *) vmalloc(vmblock, VEKTOR, n+1, 0);
if (! vmcomplete(vmblock))
**
LogError ("No Memory", 0, __FILE__, __LINE__);
return 1; //not enough memory
}
for (i=0; i<=n; i++)
for (j=0; j<=n; j++)
A[i][j] = 0.0;
// define matrix A column by column
A[1][1]=1.0; A[2][1]=-1.0;
A[1][2]=0.5; A[2][2]=0.5; A[3][2]=-1.0;
A[2][3]=1.0; A[4][3]=-1.0;
A[3][4]=2.5; A[4][4]=-1.5; A[5][4]=-1.0;
A[4][5]=4.0; A[5][5]=-3.0; A[6][5]=-1.0;
A[5][6]=5.5; A[6][6]=-4.5; A[7][6]=-1.0;
A[6][7]=7.0; A[7][7]=-6.0; A[8][7]=-1.0;
A[7][8]=8.5; A[8][8]=-7.5; A[9][8]=-1.0;
A[8][9]=10.0; A[9][9]=-9.0;
//call LU decomposition
rc = LUDCMP(A,n,INDX,&id);
//calculate determinant and display
det = id;
if (rc==0) **
for (i=1; i<=n; i++) det *= A[i][i];
printf("\n Determinant = %f\n\n", det);
return 0;
}
else **
printf("\n Error in LU decomposition.\n\n");
return 2;
}
} // main
// end of file deter1.cpp