Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
This community is designed for sharing of public information. Please do not share Intel or third-party confidential information here.
6588 Discussions

## Matrix Determinant Beginner
322 Views

Hi,

there is any way to calculate the determinant of a matrix in CSR3AC or COO format?

Thanks.

15 Replies Black Belt
322 Views
Is finding the determinant the only objective, or is that just a side result in some other computation? If the latter, what does that computation consist of? What is the typical size of the matrix? Beginner
322 Views
No, the objective is to obtain the inverse matrix.

The size of the matrix can go from 4x4 to 2000x2000 elements, where about the 20% of the values is diferent from zero.

The matrix will always be square. Employee
322 Views
Hi,
So you want to getdeterminant of sparse matrix or inverse of it? If you want to get inverse matrix you could obtain it by PARDISO functionality with many rhs.
With best regards, Beginner
322 Views
I want to get the determinant of a sparse matrix, but I only have found the functions:

mkl_dcsrsm
mkl_dcoosm

where I can assign alpha = 1 and B = to theidentity matrix in order to get the inverse of the matrix A. But the documentation of these functions does not say anything about what happens when the determinant of A is zero. That is the reason I was looking for a function to get the determinant.

I did not know about PARDISO to obtain the inverse. Is this the best option for the type of matrices I am working with?

Thanks Black Belt
322 Views
You do not need to compute the determinant to establish whether the matrix is singular. The routines used for matrix factorization will return with a suitable error code when singularity (to a user specified tolerance) is detected.

Please state what you propose to do with the inverse. Quite often, if the solution of sparse systems of linear equations with multiple right hand sides is needed, that can be obtained without using the inverse. In fact, there are many reasons why the inverse should not be explicitly formed, except in a small number of special cases. Employee
322 Views
Hi,
Please look this topic: (http://software.intel.com/en-us/forums/showthread.php?t=80373), there is a way to compute invert sparse matrix. Using dss_statistics routine (lookIntel Math Kernel Library ReferenceManualfor more details)you could obtain determinant of initial matrix.
With best regards, Beginner
322 Views

I need the inverse for the next operation:

M = -inv(A)*B

where A is of size NxN and B is of size NxM. A and B are sparse matrices with less than 20% of elements diferent from zero.

A is a submatrix from a sparse matrix C of size NxK (K>N). WhenA is singular, I have to obtain another submatrix from C, and try again.

Thanks Black Belt
322 Views
Fair enough. The singularity of A, if true, will manifest itself during the factorization of A. Note that B is not involved in this question at all. The inverse of A should not be sought. Here is a sketch of the algorithm:

1. Choose a submatrix of C as A.
2. Attempt the sparse factorization of A.
3. If Step-2 fails,
3a. Go to Step-1, choosing a different submatrix.
else
3b. Solve A X = B using the factors of A formed in Step-2.
-- done -- Beginner
322 Views

I am testing with pardiso algorithm.

I execute the phases 11, 22, and 33. When the matrix A is singular,
the phase 33 returns the values IND in the array X.

But, how can I know in the phase 22 if the matrix A is singular? Employee
322 Views
Hi,
The simplest way is to use dss interface instead of PARDISO. Then, as I mentioned before, dss_statistics routine could return determinant of initial matrix. To use PARDISO interface i need to know the value of mtype in you program, could you provide it?
With best regards, Beginner
322 Views
Hi Alexander,

the value of mtype is 11.

Thanks. Employee
322 Views
Hi,
it is not clear for me whay you got IND values in solution vector. Could you send small testcase and linkline to me to reproduce this problem on my side? By the way, in your term, as I understood you want to find matrix C with nonzero determinant as submatrix A. Could it happened that there is not such submatrix? (rang(A) < min (N,M))
With best regards, Beginner
322 Views
Hi,

this is a code where I obtain the values IND in the solution vector:

```[cpp]#include
#include

using std::cout;
using std::endl;

int main()
{

int *pt = new int;
for(int i=0; i<64; i++)	pt = 0;

//Maximal number of factors with identical nonzero sparsity structure
int maxfct = 1;

//Actual matrix for the solution phase: 1 <= mnum <= maxfct
int mnum = 1;

//Defines the matrix type
//11: real and unsymmetric matrix
int mtype = 11;

//Controls the execution of the solver
int phase = 11;

//Number of equations
int n = 2;

//Non-zero values of the coefficient matrix A

//3  1	= A
//0  0
double *a = new double;
a = 3.0;
a = 1.0;
a = 0.0;

//Row Index
int *ia = new int[n+1];
ia = 1;
ia = 3;
ia = 4;

//Columns
int *ja = new int;
ja = 1;
ja = 2;
ja = 1;

//Permutation vector of size n
int *perm = new int;
for(int i=0; i = 0;

//Number of right-hand sides that need to be solved for
int nrhs = 2;

//Pass parameters to Pardiso
int iparm;
for(int i=0; i<64; i++)	iparm = 0;

iparm = 0;		// Solver default if iparm = 0
iparm = 2;		// Fill-in reordering from METIS
// Numbers of processors, value of OMP_NUM_THREADS
iparm = 1;		// Not in use
iparm = 0;		// No iterative-direct algorithm
iparm = 0;		// No user fill-in reducing permutation
iparm = 0;		// 0: Write solution into x / 1: The solver stores the solution on the righthand side b
iparm = 0;		// Not in use
iparm = 2;		// Max numbers of iterative refinement steps
iparm = 0;		// Not in use
iparm = 13;		// Perturb the pivot elements with 1E-13
iparm = 1;		// Use nonsymmetric permutation and scaling MPS
iparm = 0;		// Not in use
iparm = 1;		// Maximum weighted matching algorithm is switched-on (default for non-symmetric)
iparm = 0;		// Output: Number of perturbed pivots
iparm = 0;		// Not in use
iparm = 0;		// Not in use
iparm = 0;		// Not in use
iparm = -1;		// Output: Number of nonzeros in the factor LU
iparm = -1;		// Output: Mflops for LU factorization
//iparm = 0;		// Output: Numbers of CG Iterations

//Message level information
//0: Pardiso generates no output / 1: Pardiso prints statistical information to the screen
int msglvl = 1;

//Right-hand side
// 5	= B
// 6
double *b = new double[n*nrhs];
for(int i=0; i = 0.0;

b = 5.0;
b = 6.0;

//Output if iparm(6)=0
double *x =  new double[n*nrhs];
for(int i=0; i = 0.0;

//Error Indicator
int error = 0;

cout << "PARDISO PHASE: " << phase << endl;

pardiso(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);

cout << "Error = " << error << endl;

cout << "b: ";
for(int i=0; i << " ";
cout << endl;

cout << "x: ";
for(int i=0; i << " ";
cout << endl;

phase = 22;

cout << "PARDISO PHASE: " << phase << endl;

pardiso(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);

cout << "Error =  " << error << endl;

cout << "b: ";
for(int i=0; i << " ";
cout << endl;

cout << "x: ";
for(int i=0; i << " ";
cout << endl;

phase = 33;

cout << "PARDISO PHASE: " << phase << endl;

pardiso(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);

cout << "Error = " << error << endl;

cout << "b: ";
for(int i=0; i << " ";
cout << endl;

cout << "x: ";
for(int i=0; i << " ";
cout << endl;

}[/cpp]```

The matrix A is singular:

3 1
00

Respect to your question, always exist a submatrix of C that is not singular.

Thanks Employee
322 Views
Hi,
I've reproduced your issue, I will try to resolve it asap.
With best regards,  