本文整理汇总了C#中matinvreport类的典型用法代码示例。如果您正苦于以下问题：C# matinvreport类的具体用法？C# matinvreport怎么用？C# matinvreport使用的例子？那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
matinvreport类属于命名空间，在下文中一共展示了matinvreport类的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C#代码示例。
示例1: A
/*************************************************************************
Inversion of a Hermitian positive definite matrix which is given
by Cholesky decomposition.
Input parameters:
A - Cholesky decomposition of the matrix to be inverted:
A=U’*U or A = L*L'.
Output of HPDMatrixCholesky subroutine.
N - size of matrix A (optional) :
* if given, only principal NxN submatrix is processed and
overwritten. other elements are unchanged.
* if not given, size is automatically determined from
matrix size (A must be square matrix)
IsUpper - storage type (optional):
* if True, symmetric matrix A is given by its upper
triangle, and the lower triangle isn’t used/changed by
function
* if False, symmetric matrix A is given by its lower
triangle, and the upper triangle isn’t used/changed by
function
* if not given, lower half is used.
Output parameters:
Info - return code, same as in RMatrixLUInverse
Rep - solver report, same as in RMatrixLUInverse
A - inverse of matrix A, same as in RMatrixLUInverse
-- ALGLIB routine --
10.02.2010
Bochkanov Sergey
*************************************************************************/
public static void hpdmatrixcholeskyinverse(ref complex[,] a,
int n,
bool isupper,
ref int info,
matinvreport rep)
{
int i = 0;
int j = 0;
matinvreport rep2 = new matinvreport();
complex[] tmp = new complex[0];
bool f = new bool();
info = 0;
ap.assert(n>0, "HPDMatrixCholeskyInverse: N<=0!");
ap.assert(ap.cols(a)>=n, "HPDMatrixCholeskyInverse: cols(A)<N!");
ap.assert(ap.rows(a)>=n, "HPDMatrixCholeskyInverse: rows(A)<N!");
f = true;
for(i=0; i<=n-1; i++)
{
f = (f & math.isfinite(a[i,i].x)) & math.isfinite(a[i,i].y);
}
ap.assert(f, "HPDMatrixCholeskyInverse: A contains infinite or NaN values!");
info = 1;
//
// calculate condition numbers
//
rep.r1 = rcond.hpdmatrixcholeskyrcond(a, n, isupper);
rep.rinf = rep.r1;
if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) | (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
{
if( isupper )
{
for(i=0; i<=n-1; i++)
{
for(j=i; j<=n-1; j++)
{
a[i,j] = 0;
}
}
}
else
{
for(i=0; i<=n-1; i++)
{
for(j=0; j<=i; j++)
{
a[i,j] = 0;
}
}
}
rep.r1 = 0;
rep.rinf = 0;
info = -3;
return;
}
//
// Inverse
//
tmp = new complex[n];
hpdmatrixcholeskyinverserec(ref a, 0, n, isupper, ref tmp);
}
开发者ID:Ring-r，项目名称:opt，代码行数:95，代码来源:linalg.cs
示例2: inverted
/*************************************************************************
Inversion of a Hermitian positive definite matrix.
Given an upper or lower triangle of a Hermitian positive definite matrix,
the algorithm generates matrix A^-1 and saves the upper or lower triangle
depending on the input.
Input parameters:
A - matrix to be inverted (upper or lower triangle).
Array with elements [0..N-1,0..N-1].
N - size of matrix A (optional) :
* if given, only principal NxN submatrix is processed and
overwritten. other elements are unchanged.
* if not given, size is automatically determined from
matrix size (A must be square matrix)
IsUpper - storage type (optional):
* if True, symmetric matrix A is given by its upper
triangle, and the lower triangle isn’t used/changed by
function
* if False, symmetric matrix A is given by its lower
triangle, and the upper triangle isn’t used/changed by
function
* if not given, both lower and upper triangles must be
filled.
Output parameters:
Info - return code, same as in RMatrixLUInverse
Rep - solver report, same as in RMatrixLUInverse
A - inverse of matrix A, same as in RMatrixLUInverse
-- ALGLIB routine --
10.02.2010
Bochkanov Sergey
*************************************************************************/
public static void hpdmatrixinverse(ref complex[,] a,
int n,
bool isupper,
ref int info,
matinvreport rep)
{
info = 0;
ap.assert(n>0, "HPDMatrixInverse: N<=0!");
ap.assert(ap.cols(a)>=n, "HPDMatrixInverse: cols(A)<N!");
ap.assert(ap.rows(a)>=n, "HPDMatrixInverse: rows(A)<N!");
ap.assert(apserv.apservisfinitectrmatrix(a, n, isupper), "HPDMatrixInverse: A contains infinite or NaN values!");
info = 1;
if( trfac.hpdmatrixcholesky(ref a, n, isupper) )
{
hpdmatrixcholeskyinverse(ref a, n, isupper, ref info, rep);
}
else
{
info = -3;
}
}
开发者ID:Ring-r，项目名称:opt，代码行数:56，代码来源:linalg.cs
示例3: matrix
/*************************************************************************
Inversion of a matrix given by its LU decomposition.
INPUT PARAMETERS:
A - LU decomposition of the matrix (output of RMatrixLU subroutine).
Pivots - table of permutations which were made during the LU decomposition
(the output of RMatrixLU subroutine).
N - size of matrix A.
OUTPUT PARAMETERS:
Info - return code:
* -3 A is singular, or VERY close to singular.
it is filled by zeros in such cases.
* -1 N<=0 was passed, or incorrect Pivots was passed
* 1 task is solved (but matrix A may be ill-conditioned,
check R1/RInf parameters for condition numbers).
Rep - solver report, see below for more info
A - inverse of matrix A.
Array whose indexes range within [0..N-1, 0..N-1].
SOLVER REPORT
Subroutine sets following fields of the Rep structure:
* R1 reciprocal of condition number: 1/cond(A), 1-norm.
* RInf reciprocal of condition number: 1/cond(A), inf-norm.
-- ALGLIB routine --
05.02.2010
Bochkanov Sergey
*************************************************************************/
public static void rmatrixluinverse(ref double[,] a,
ref int[] pivots,
int n,
ref int info,
ref matinvreport rep)
{
double[] work = new double[0];
int i = 0;
int j = 0;
int k = 0;
double v = 0;
info = 1;
//
// Quick return if possible
//
if (n == 0)
{
info = -1;
return;
}
for (i = 0; i <= n - 1; i++)
{
if (pivots[i] > n - 1 | pivots[i] < i)
{
info = -1;
return;
}
}
//
// calculate condition numbers
//
rep.r1 = rcond.rmatrixlurcond1(ref a, n);
rep.rinf = rcond.rmatrixlurcondinf(ref a, n);
if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold()))
{
for (i = 0; i <= n - 1; i++)
{
for (j = 0; j <= n - 1; j++)
{
a[i, j] = 0;
}
}
rep.r1 = 0;
rep.rinf = 0;
info = -3;
return;
}
//
// Call cache-oblivious code
//
work = new double[n];
rmatrixluinverserec(ref a, 0, n, ref work, ref info, ref rep);
//
// apply permutations
//
for (i = 0; i <= n - 1; i++)
{
for (j = n - 2; j >= 0; j--)
{
k = pivots[j];
v = a[i, j];
a[i, j] = a[i, k];
a[i, k] = v;
}
}
//.........这里部分代码省略.........
开发者ID:mortenbakkedal，项目名称:SharpMath，代码行数:101，代码来源:matinv.cs
示例4: _pexec_rmatrixinverse
/*************************************************************************
Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
*************************************************************************/
public static void _pexec_rmatrixinverse(ref double[,] a,
int n,
ref int info,
matinvreport rep)
{
rmatrixinverse(ref a,n,ref info,rep);
}
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:10，代码来源:linalg.cs
示例5: cmatrixinverse
/*************************************************************************
Inversion of a general matrix.
Input parameters:
A - matrix, array[0..N-1,0..N-1].
N - size of A.
Output parameters:
Info - return code, same as in RMatrixLUInverse
Rep - solver report, same as in RMatrixLUInverse
A - inverse of matrix A, same as in RMatrixLUInverse
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*************************************************************************/
public static void cmatrixinverse(ref AP.Complex[,] a,
int n,
ref int info,
ref matinvreport rep)
{
int[] pivots = new int[0];
trfac.cmatrixlu(ref a, n, n, ref pivots);
cmatrixluinverse(ref a, ref pivots, n, ref info, ref rep);
}
开发者ID:mortenbakkedal，项目名称:SharpMath，代码行数:25，代码来源:matinv.cs
示例6: inverted
/*************************************************************************
Inversion of a Hermitian positive definite matrix.
Given an upper or lower triangle of a Hermitian positive definite matrix,
the algorithm generates matrix A^-1 and saves the upper or lower triangle
depending on the input.
Input parameters:
A - matrix to be inverted (upper or lower triangle).
Array with elements [0..N-1,0..N-1].
N - size of matrix A.
IsUpper - storage format.
If IsUpper = True, then the upper triangle of matrix A is
given, otherwise the lower triangle is given.
Output parameters:
Info - return code, same as in RMatrixLUInverse
Rep - solver report, same as in RMatrixLUInverse
A - inverse of matrix A, same as in RMatrixLUInverse
-- ALGLIB routine --
10.02.2010
Bochkanov Sergey
*************************************************************************/
public static void hpdmatrixinverse(ref AP.Complex[,] a,
int n,
bool isupper,
ref int info,
ref matinvreport rep)
{
if (n < 1)
{
info = -1;
return;
}
info = 1;
if (trfac.hpdmatrixcholesky(ref a, n, isupper))
{
hpdmatrixcholeskyinverse(ref a, n, isupper, ref info, ref rep);
}
else
{
info = -3;
}
}
开发者ID:mortenbakkedal，项目名称:SharpMath，代码行数:45，代码来源:matinv.cs
示例7: _pexec_cmatrixtrinverse
/*************************************************************************
Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
*************************************************************************/
public static void _pexec_cmatrixtrinverse(ref complex[,] a,
int n,
bool isupper,
bool isunit,
ref int info,
matinvreport rep)
{
cmatrixtrinverse(ref a,n,isupper,isunit,ref info,rep);
}
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:12，代码来源:linalg.cs
示例8: hpdmatrixcholeskyinverserec
/*************************************************************************
Recursive subroutine for HPD inversion.
-- ALGLIB routine --
10.02.2010
Bochkanov Sergey
*************************************************************************/
private static void hpdmatrixcholeskyinverserec(ref AP.Complex[,] a,
int offs,
int n,
bool isupper,
ref AP.Complex[] tmp)
{
int i = 0;
int j = 0;
AP.Complex v = 0;
int n1 = 0;
int n2 = 0;
int info2 = 0;
matinvreport rep2 = new matinvreport();
int i_ = 0;
int i1_ = 0;
if (n < 1)
{
return;
}
//
// Base case
//
if (n <= ablas.ablascomplexblocksize(ref a))
{
cmatrixtrinverserec(ref a, offs, n, isupper, false, ref tmp, ref info2, ref rep2);
if (isupper)
{
//
// Compute the product U * U'.
// NOTE: we never assume that diagonal of U is real
//
for (i = 0; i <= n - 1; i++)
{
if (i == 0)
{
//
// 1x1 matrix
//
a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
}
else
{
//
// (I+1)x(I+1) matrix,
//
// ( A11 A12 ) ( A11^H ) ( A11*A11^H+A12*A12^H A12*A22^H )
// ( ) * ( ) = ( )
// ( A22 ) ( A12^H A22^H ) ( A22*A12^H A22*A22^H )
//
// A11 is IxI, A22 is 1x1.
//
i1_ = (offs) - (0);
for (i_ = 0; i_ <= i - 1; i_++)
{
tmp[i_] = AP.Math.Conj(a[i_ + i1_, offs + i]);
}
for (j = 0; j <= i - 1; j++)
{
v = a[offs + j, offs + i];
i1_ = (j) - (offs + j);
for (i_ = offs + j; i_ <= offs + i - 1; i_++)
{
a[offs + j, i_] = a[offs + j, i_] + v * tmp[i_ + i1_];
}
}
v = AP.Math.Conj(a[offs + i, offs + i]);
for (i_ = offs; i_ <= offs + i - 1; i_++)
{
a[i_, offs + i] = v * a[i_, offs + i];
}
a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
}
}
}
else
{
//
// Compute the product L' * L
// NOTE: we never assume that diagonal of L is real
//
for (i = 0; i <= n - 1; i++)
{
if (i == 0)
{
//
// 1x1 matrix
//.........这里部分代码省略.........
开发者ID:mortenbakkedal，项目名称:SharpMath，代码行数:101，代码来源:matinv.cs
示例9: inverse
/*************************************************************************
Triangular matrix inverse (real)
The subroutine inverts the following types of matrices:
* upper triangular
* upper triangular with unit diagonal
* lower triangular
* lower triangular with unit diagonal
In case of an upper (lower) triangular matrix, the inverse matrix will
also be upper (lower) triangular, and after the end of the algorithm, the
inverse matrix replaces the source matrix. The elements below (above) 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.
COMMERCIAL EDITION OF ALGLIB:
! Commercial version of ALGLIB includes two important improvements of
! this function, which can be used from C++ and C#:
! * Intel MKL support (lightweight Intel MKL is shipped with ALGLIB)
! * multicore support
!
! Intel MKL gives approximately constant (with respect to number of
! worker threads) acceleration factor which depends on CPU being used,
! problem size and "baseline" ALGLIB edition which is used for
! comparison.
!
! Say, on SSE2-capable CPU with N=1024, HPC ALGLIB will be:
! * about 2-3x faster than ALGLIB for C++ without MKL
! * about 7-10x faster than "pure C#" edition of ALGLIB
! Difference in performance will be more striking on newer CPU's with
! support for newer SIMD instructions. Generally, MKL accelerates any
! problem whose size is at least 128, with best efficiency achieved for
! N's larger than 512.
!
! Commercial edition of ALGLIB also supports multithreaded acceleration
! of this function. We should note that triangular inverse is harder to
! parallelize than, say, matrix-matrix product - this algorithm has
! many internal synchronization points which can not be avoided. However
! parallelism starts to be profitable starting from N=1024, achieving
! near-linear speedup for N=4096 or higher.
!
! In order to use multicore features you have to:
! * use commercial version of ALGLIB
! * call this function with "smp_" prefix, which indicates that
! multicore code will be used (for multicore support)
!
! We recommend you to read 'Working with commercial version' section of
! ALGLIB Reference Manual in order to find out how to use performance-
! related features provided by commercial edition of ALGLIB.
Input parameters:
A - matrix, array[0..N-1, 0..N-1].
N - size of matrix A (optional) :
* if given, only principal NxN submatrix is processed and
overwritten. other elements are unchanged.
* if not given, size is automatically determined from
matrix size (A must be square matrix)
IsUpper - True, if the matrix is upper triangular.
IsUnit - diagonal type (optional):
* if True, matrix has unit diagonal (a[i,i] are NOT used)
* if False, matrix diagonal is arbitrary
* if not given, False is assumed
Output parameters:
Info - same as for RMatrixLUInverse
Rep - same as for RMatrixLUInverse
A - same as for RMatrixLUInverse.
-- ALGLIB --
Copyright 05.02.2010 by Bochkanov Sergey
*************************************************************************/
public static void rmatrixtrinverse(ref double[,] a,
int n,
bool isupper,
bool isunit,
ref int info,
matinvreport rep)
{
int i = 0;
int j = 0;
double[] tmp = new double[0];
apserv.sinteger sinfo = new apserv.sinteger();
info = 0;
alglib.ap.assert(n>0, "RMatrixTRInverse: N<=0!");
alglib.ap.assert(alglib.ap.cols(a)>=n, "RMatrixTRInverse: cols(A)<N!");
alglib.ap.assert(alglib.ap.rows(a)>=n, "RMatrixTRInverse: rows(A)<N!");
alglib.ap.assert(apserv.isfinitertrmatrix(a, n, isupper), "RMatrixTRInverse: A contains infinite or NaN values!");
//
// calculate condition numbers
//
rep.r1 = rcond.rmatrixtrrcond1(a, n, isupper, isunit);
rep.rinf = rcond.rmatrixtrrcondinf(a, n, isupper, isunit);
if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) || (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
{
//.........这里部分代码省略.........
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:101，代码来源:linalg.cs
示例10: _pexec_rmatrixtrinverse
/*************************************************************************
Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
*************************************************************************/
public static void _pexec_rmatrixtrinverse(ref double[,] a,
int n,
bool isupper,
bool isunit,
ref int info,
matinvreport rep)
{
rmatrixtrinverse(ref a,n,isupper,isunit,ref info,rep);
}
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:12，代码来源:linalg.cs
示例11: _pexec_hpdmatrixinverse
/*************************************************************************
Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
*************************************************************************/
public static void _pexec_hpdmatrixinverse(ref complex[,] a,
int n,
bool isupper,
ref int info,
matinvreport rep)
{
hpdmatrixinverse(ref a,n,isupper,ref info,rep);
}
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:11，代码来源:linalg.cs
示例12: _pexec_spdmatrixinverse
/*************************************************************************
Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
*************************************************************************/
public static void _pexec_spdmatrixinverse(ref double[,] a,
int n,
bool isupper,
ref int info,
matinvreport rep)
{
spdmatrixinverse(ref a,n,isupper,ref info,rep);
}
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:11，代码来源:linalg.cs
示例13: _pexec_cmatrixinverse
/*************************************************************************
Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
*************************************************************************/
public static void _pexec_cmatrixinverse(ref complex[,] a,
int n,
ref int info,
matinvreport rep)
{
cmatrixinverse(ref a,n,ref info,rep);
}
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:10，代码来源:linalg.cs
示例14: inverse
/*************************************************************************
Triangular matrix inverse (complex)
The subroutine inverts the following types of matrices:
* upper triangular
* upper triangular with unit diagonal
* lower triangular
* lower triangular with unit diagonal
In case of an upper (lower) triangular matrix, the inverse matrix will
also be upper (lower) triangular, and after the end of the algorithm, the
inverse matrix replaces the source matrix. The elements below (above) 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:
A - matrix, array[0..N-1, 0..N-1].
N - size of matrix A (optional) :
* if given, only principal NxN submatrix is processed and
overwritten. other elements are unchanged.
* if not given, size is automatically determined from
matrix size (A must be square matrix)
IsUpper - True, if the matrix is upper triangular.
IsUnit - diagonal type (optional):
* if True, matrix has unit diagonal (a[i,i] are NOT used)
* if False, matrix diagonal is arbitrary
* if not given, False is assumed
Output parameters:
Info - same as for RMatrixLUInverse
Rep - same as for RMatrixLUInverse
A - same as for RMatrixLUInverse.
-- ALGLIB --
Copyright 05.02.2010 by Bochkanov Sergey
*************************************************************************/
public static void cmatrixtrinverse(ref complex[,] a,
int n,
bool isupper,
bool isunit,
ref int info,
matinvreport rep)
{
int i = 0;
int j = 0;
complex[] tmp = new complex[0];
info = 0;
ap.assert(n>0, "CMatrixTRInverse: N<=0!");
ap.assert(ap.cols(a)>=n, "CMatrixTRInverse: cols(A)<N!");
ap.assert(ap.rows(a)>=n, "CMatrixTRInverse: rows(A)<N!");
ap.assert(apserv.apservisfinitectrmatrix(a, n, isupper), "CMatrixTRInverse: A contains infinite or NaN values!");
info = 1;
//
// calculate condition numbers
//
rep.r1 = rcond.cmatrixtrrcond1(a, n, isupper, isunit);
rep.rinf = rcond.cmatrixtrrcondinf(a, n, isupper, isunit);
if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) | (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
{
for(i=0; i<=n-1; i++)
{
for(j=0; j<=n-1; j++)
{
a[i,j] = 0;
}
}
rep.r1 = 0;
rep.rinf = 0;
info = -3;
return;
}
//
// Invert
//
tmp = new complex[n];
cmatrixtrinverserec(ref a, 0, n, isupper, isunit, ref tmp, ref info, rep);
}
开发者ID:Ring-r，项目名称:opt，代码行数:83，代码来源:linalg.cs
示例15: spdmatrixcholeskyinverserec
/*************************************************************************
Recursive subroutine for SPD inversion.
NOTE: this function expects that matris is strictly positive-definite.
-- ALGLIB routine --
10.02.2010
Bochkanov Sergey
*************************************************************************/
public static void spdmatrixcholeskyinverserec(ref double[,] a,
int offs,
int n,
bool isupper,
ref double[] tmp)
{
int i = 0;
int j = 0;
double v = 0;
int n1 = 0;
int n2 = 0;
apserv.sinteger sinfo2 = new apserv.sinteger();
matinvreport rep2 = new matinvreport();
int i_ = 0;
int i1_ = 0;
if( n<1 )
{
return;
}
//
// Base case
//
if( n<=ablas.ablasblocksize(a) )
{
sinfo2.val = 1;
rmatrixtrinverserec(a, offs, n, isupper, false, tmp, sinfo2, rep2);
alglib.ap.assert(sinfo2.val>0, "SPDMatrixCholeskyInverseRec: integrity check failed");
if( isupper )
{
//
// Compute the product U * U'.
// NOTE: we never assume that diagonal of U is real
//
for(i=0; i<=n-1; i++)
{
if( i==0 )
{
//
// 1x1 matrix
//
a[offs+i,offs+i] = math.sqr(a[offs+i,offs+i]);
}
else
{
//
// (I+1)x(I+1) matrix,
//
// ( A11 A12 ) ( A11^H ) ( A11*A11^H+A12*A12^H A12*A22^H )
// ( ) * ( ) = ( )
// ( A22 ) ( A12^H A22^H ) ( A22*A12^H A22*A22^H )
//
// A11 is IxI, A22 is 1x1.
//
i1_ = (offs) - (0);
for(i_=0; i_<=i-1;i_++)
{
tmp[i_] = a[i_+i1_,offs+i];
}
for(j=0; j<=i-1; j++)
{
v = a[offs+j,offs+i];
i1_ = (j) - (offs+j);
for(i_=offs+j; i_<=offs+i-1;i_++)
{
a[offs+j,i_] = a[offs+j,i_] + v*tmp[i_+i1_];
}
}
v = a[offs+i,offs+i];
for(i_=offs; i_<=offs+i-1;i_++)
{
a[i_,offs+i] = v*a[i_,offs+i];
}
a[offs+i,offs+i] = math.sqr(a[offs+i,offs+i]);
}
}
}
else
{
//
// Compute the product L' * L
// NOTE: we never assume that diagonal of L is real
//
for(i=0; i<=n-1; i++)
{
if( i==0 )
//.........这里部分代码省略.........
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:101，代码来源:linalg.cs
示例16: cmatrixluinverserec
private static void cmatrixluinverserec(ref AP.Complex[,] a,
int offs,
int n,
ref AP.Complex[] work,
ref int info,
ref matinvreport rep)
{
int i = 0;
int iws = 0;
int j = 0;
int jb = 0;
int jj = 0;
int jp = 0;
int k = 0;
AP.Complex v = 0;
int n1 = 0;
int n2 = 0;
int i_ = 0;
int i1_ = 0;
if (n < 1)
{
info = -1;
return;
}
//
// Base case
//
if (n <= ablas.ablascomplexblocksize(ref a))
{
//
// Form inv(U)
//
cmatrixtrinverserec(ref a, offs, n, true, false, ref work, ref info, ref rep);
if (info <= 0)
{
return;
}
//
// Solve the equation inv(A)*L = inv(U) for inv(A).
//
for (j = n - 1; j >= 0; j--)
{
//
// Copy current column of L to WORK and replace with zeros.
//
for (i = j + 1; i <= n - 1; i++)
{
work[i] = a[offs + i, offs + j];
a[offs + i, offs + j] = 0;
}
//
// Compute current column of inv(A).
//
if (j < n - 1)
{
for (i = 0; i <= n - 1; i++)
{
i1_ = (j + 1) - (offs + j + 1);
v = 0.0;
for (i_ = offs + j + 1; i_ <= offs + n - 1; i_++)
{
v += a[offs + i, i_] * work[i_ + i1_];
}
a[offs + i, offs + j] = a[offs + i, offs + j] - v;
}
}
}
return;
}
//
// Recursive code:
//
// ( L1 ) ( U1 U12 )
// A = ( ) * ( )
// ( L12 L2 ) ( U2 )
//
// ( W X )
// A^-1 = ( )
// ( Y Z )
//
ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
System.Diagnostics.Debug.Assert(n2 > 0, "LUInverseRec: internal error!");
//
// X := inv(U1)*U12*inv(U2)
//
ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, true, false, 0, ref a, offs, offs + n1);
ablas.cmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, true, false, 0, ref a, offs, offs + n1);
//
// Y := inv(L2)*L12*inv(L1)
//
ablas.cmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, false, true, 0, ref a, offs + n1, offs);
//.........这里部分代码省略.........
开发者ID:mortenbakkedal，项目名称:SharpMath，代码行数:101，代码来源:matinv.cs
示例17: rmatrixtrinverserec
/*************************************************************************
Triangular matrix inversion, recursive subroutine
NOTE: this function sets Info on failure, leaves it unchanged on success.
NOTE: only Tmp[Offs:Offs+N-1] is modified, other entries of the temporary array are not modified
-- ALGLIB --
05.02.2010, Bochkanov Sergey.
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
February 29, 1992.
*************************************************************************/
private static void rmatrixtrinverserec(double[,] a,
int offs,
int n,
bool isupper,
bool isunit,
double[] tmp,
apserv.sinteger info,
matinvreport rep)
{
int n1 = 0;
int n2 = 0;
int i = 0;
int j = 0;
double v = 0;
double ajj = 0;
int i_ = 0;
if( n<1 )
{
info.val = -1;
return;
}
//
// Base case
//
if( n<=ablas.ablasblocksize(a) )
{
if( isupper )
{
//
// Compute inverse of upper triangular matrix.
//
for(j=0; j<=n-1; j++)
{
if( !isunit )
{
if( (double)(a[offs+j,offs+j])==(double)(0) )
{
info.val = -3;
return;
}
a[offs+j,offs+j] = 1/a[offs+j,offs+j];
ajj = -a[offs+j,offs+j];
}
else
{
ajj = -1;
}
//
// Compute elements 1:j-1 of j-th column.
//
if( j>0 )
{
for(i_=offs+0; i_<=offs+j-1;i_++)
{
tmp[i_] = a[i_,offs+j];
}
for(i=0; i<=j-1; i++)
{
if( i<j-1 )
{
v = 0.0;
for(i_=offs+i+1; i_<=offs+j-1;i_++)
{
v += a[offs+i,i_]*tmp[i_];
}
}
else
{
v = 0;
}
if( !isunit )
{
a[offs+i,offs+j] = v+a[offs+i,offs+i]*tmp[offs+i];
}
else
{
a[offs+i,offs+j] = v+tmp[offs+i];
}
}
for(i_=offs+0; i_<=offs+j-1;i_++)
{
a[i_,offs+j] = ajj*a[i_,offs+j];
}
//.........这里部分代码省略.........
开发者ID:Kerbas-ad-astra，项目名称:MechJeb2，代码行数:101，代码来源:linalg.cs
示例18: rmatrixinverse
/*************************************************************************
Inversion of a general matrix.
Input parameters:
A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
N - size of matrix A.
Output parameters:
Info - return code, same as in RMatrixLUInverse
Rep - solver report, same as in RMatrixLUInverse
A - inverse of matrix A, same as in RMatrixLUInverse
Result:
True, if the matrix is not singular.
False, if the matrix is singular.
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*********************************************************************
|
客服电话
APP下载
官方微信
请发表评论