mkl_ddnscsr problems

[Jared_W_]

Hello everyone

I am trying to use the ddnscsr to build the arrays a, ia, ja which will then be used for the pardiso or other routines which need CSR formatting. 

Basically to make sure I understand calling the routine I have been trying to implement it on the nonsymmetric matrix listed in the sparse storage format page, but it isnt working and I cannot figure out why.

    //
    //                 |   1    -1    0   -3   0   |
    //                 |  -2     5    0    0   0   |
    //   A    =           |   0     0    4    6   4   |
    //                 |  -4     0    2    7   0   |
    //                 |   0     8    0    0   -5  |
    //
    //  

    vector<vector<double>> A(5, vector<double>(5, 0));

    A[0][0] = 1.0;
    A[0][1] = -1.0;
    A[0][2] = 0.0;
    A[0][3] = -3.0;
    A[0][4] = 0.0;

    A[1][0] = -2.0;
    A[1][1] = 5.0;
    A[1][2] = 0.0;
    A[1][3] = 0.0;
    A[1][4] = 0.0;

    A[2][0] = 0.0;
    A[2][1] = 0.0;
    A[2][2] = 4.0;
    A[2][3] = 6.0;
    A[2][4] = 4.0;

    A[3][0] = -4.0;
    A[3][1] = 0.0;
    A[3][2] = 2.0;
    A[3][3] = 7.0;
    A[3][4] = 0.0;

    A[4][0] = 0.0;
    A[4][1] = 8.0;
    A[4][2] = 0.0;
    A[4][3] = 0.0;
    A[4][4] = -5.0;

    // ia[5] = {0, 3, 5, 8, 11, 14}
    // ja[13] = {0, 1, 3, 0, 1, 2, 3, 4, 0, 2, 3, 1, 4}
    // a[13] = {1, -1, -3, -2, 5, 4, 6, 4, -4, 2, 7, 8, -5}

            double *solution, *adns, *a, *b, *x;
            int rows, cols, c, i, j, k;
            MKL_INT m, n, lda, info;
            lda = 1;
            m = 5;
            n = m;
            int nzmax = 13;

            adns = (double*)mkl_malloc(nzmax*sizeof(double),16);
            if (adns == NULL) 
            {
                cout << ">>> error allocating adns" << endl;
                return (0);
            }
            a = (double*)mkl_malloc(nzmax*sizeof(double),16);
            if (a == NULL) 
            {
                cout << ">>> error allocating a" << endl;
                return (0);
            }

            MKL_INT *ia, *ja;
            ia = (MKL_INT *)mkl_malloc((m + 1)*sizeof(MKL_INT),16);
            if (ia == NULL) 
            {
                cout << ">>> error allocating a" << endl;
                return (0);
            }
            ja = (MKL_INT*)mkl_malloc(nzmax*sizeof(MKL_INT),16);
            if (ja == NULL) 
            {
                cout << ">>> error allocating a" << endl;
                return (0);
            }

            /*    Building input parameters for pardiso()    */
            /*    adns - array */
            c = 0;
            for (rows=0; rows<m; rows++) 
            {
                for (cols=0; cols<n; cols++) 
                {
                    if(A[rows][cols] != 0)
                    {
                        adns = A[rows][cols];
                        c++;
                    }
                }
            }

            MKL_INT *job;
            job = (MKL_INT*)mkl_malloc((m + 1)*sizeof(MKL_INT),16);
            if (ja == NULL) 
            {
                cout << ">>> error allocating a" << endl;
                return (0);
            }
            job[0] = 0;
            job[1] = 0;
            job[2] = 0;
            job[3] = 2;
            job[4] = 13;
            job[5] = 3;

            mkl_ddnscsr(job, &m, &n, adns, &lda, a, ja, ia, &info);

Here is my whole code I put together. the "A" matrix is the simple 5 X 5 matrix as shown above, and I even commented out the correct format as show, but for some reason it will not print out the a, ja, and ia arrays properly. Currently, the info = 2 and I cannot see why it is failing on the second row... please someone shed some light.

I have also attached the current output from this code.

[mecej4]

You have variables adns and a interchanged. The adjectives "dense" and "sparse" refer to the storage scheme, not the matrices! The "dense" array is an array of 25 elements of the matrix A, including even the zero elements. The output "sparse" matrix values in array a do not contain any zeros.

Here is a short C program to illustrate correct usage and conversion (unfortunately, mangled by the IDF formatter:by adding blank lines arbitrarily, and not displaying the backslash character).

[cpp]

#include <stdio.h>
#include <mkl.h>
int main(){
double A[5][5]={{1,-1,0,-3,0},{-2,5,0,0,0},{0,0,4,6,4},{-4,0,2,7,0},
                {0,8,0,0,-5}};
int job[6]={0,0,0,2,13,3};
int m=5,n=5,lda=5,nnz=13,i,info;
int ja[13],ia[6];
double aval[13];

mkl_ddnscsr(job,&m,&n,(double *)A,&lda,aval,ja,ia,&info);
printf("info = %d\n",info);
printf("Row pointers\n");
for(i=0; i<m+1; i++)printf("%2d %2d ",i,ia);
printf("\n Columns and values\n");
for(i=0; i<nnz; i++)printf("%2d %2d %6.1f\n",i,ja,aval);
}

[/cpp]

[Zhang_Z_Intel]

As mecej4 pointed out, adns is an input dense matrix of full storage format, its size should be m*n=25.

Another problem in your code, lda should be 5, not 1. This is the leading dimension of the matrix, which is the number of elements in the major dimension.

[Jared_W_]

Thank you guys so much, it is amazing how fast you guys respond. I'm not sure if I am not understanding something properly in the reference manual, but it mentions that adns is, "Array containing non-zero elements of the matrix A." Am I off or is the reference manual information mixed up?

link to reference manual:

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm

Zhang thanks for the quick reply and the sharp eye on my code.

Mecej4 you have been extremely helpful for me every single time I have posted I really appreciate the help.

Happy Holidays!

[Jared_W_]

One last thing and I think I will be good. I have corrected my stuff and now things are printing out but I noticed the ia array seems off. The last entry (I believe) is supposed to be the number of non zeros + 1? Which in the example from myself and mecej4 should be 14, but is printing out to be 13...Does this seem wrong or am I overlooking something obvious?

[mecej4]

The entries in the ia array are indices. Depending on the indexing convention selected, the values in that array range from 0 : n_nz or from 1 : n_nz + 1. What values did you specify in the job array?

[Jared_W_]

I specified the same job array as you. I get the same values and I feel like they are correct I was just thinking the last element of the ia array had to be (nz + 1). So for this problem nz = 13 so the last element of ia would be 14, but I'm thinking I was wrong. 

[mecej4]

 I was just thinking the last element of the ia array had to be (nz + 1)

Much of the BLAS/Lapack documentation has a Fortran heritage, and you may find vestiges of that heritage in the documentation. In older versions of Fortran programs array indices started invariably with 1 rather than zero, and some of the statements in the manuals may have that bias. On the other hand, including a parenthetical correction for 0-based indexing in the manuals would make them harder to read.