Community
cancel
Showing results for
Did you mean:
Beginner
166 Views

## Different pivots ( intel mkl vs gnu lapack) from dgesv !

i have two codes for solving linear equations, one from intel:

```
/*******************************************************************************
*  The information and material ("Material") provided below is owned by Intel
*  Corporation or its suppliers or licensors, and title to such Material remains
*  with Intel Corporation or its suppliers or licensors. The Material contains
*  proprietary information of Intel or its suppliers and licensors. The Material
*  is protected by worldwide copyright laws and treaty provisions. No part of
*  the Material may be copied, reproduced, published, uploaded, posted,
*  transmitted, or distributed in any way without Intel's prior express written
*  property rights in the Material is granted to or conferred upon you, either
*  expressly, by implication, inducement, estoppel or otherwise. Any license
*  under such intellectual property rights must be express and approved by Intel
*  in writing.
*
********************************************************************************
*/
/*
LAPACKE_dgesv Example.
======================

The program computes the solution to the system of linear
equations with a square matrix A and multiple
right-hand sides B, where A is the coefficient matrix:

6.80  -6.05  -0.45   8.32  -9.67
-2.11  -3.30   2.58   2.71  -5.14
5.66   5.36  -2.70   4.35  -7.26
5.97  -4.44   0.27  -7.17   6.08
8.23   1.08   9.04   2.14  -6.87

and B is the right-hand side matrix:

4.02  -1.56   9.81
6.19   4.00  -4.09
-8.22  -8.67  -4.57
-7.57   1.75  -8.61
-3.03   2.86   8.99

Description.
============

The routine solves for X the system of linear equations A*X = B,
where A is an n-by-n matrix, the columns of matrix B are individual
right-hand sides, and the columns of X are the corresponding
solutions.

The LU decomposition with partial pivoting and row interchanges is
used to factor A as A = P*L*U, where P is a permutation matrix, L
is unit lower triangular, and U is upper triangular. The factored
form of A is then used to solve the system of equations A*X = B.

Example Program Results.
========================

LAPACKE_dgesv (row-major, high-level) Example Program Results

Solution
-0.80  -0.39   0.96
-0.70  -0.55   0.22
0.59   0.84   1.90
1.32  -0.10   5.36
0.57   0.11   4.04

Details of LU factorization
8.23   1.08   9.04   2.14  -6.87
0.83  -6.94  -7.92   6.55  -3.99
0.69  -0.67 -14.18   7.24  -5.19
0.73   0.75   0.02 -13.82  14.19
-0.26   0.44  -0.59  -0.34  -3.43

Pivot indices
5      5      3      4      5
*/
#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"

/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda );
extern void print_int_vector( char* desc, MKL_INT n, MKL_INT* a );

/* Parameters */
#define N 3
#define NRHS 1
#define LDA N
#define LDB NRHS

/* Main program */
int main() {
/* Locals */
MKL_INT n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info;
/* Local arrays */
MKL_INT ipiv;
double a[LDA*N] = {
1,1,1,
1,1,3,
2,1,1
};
double b[LDB*N] = {
1,
2,
3
};
/* Executable statements */
printf( "LAPACKE_dgesv (row-major, high-level) Example Program Results\n" );
/* Solve the equations A*X = B */
info = LAPACKE_dgesv( LAPACK_ROW_MAJOR, n, nrhs, a, lda, ipiv,
b, ldb );
/* Check for the exact singularity */
if( info > 0 ) {
printf( "The diagonal element of the triangular factor of A,\n" );
printf( "U(%i,%i) is zero, so that A is singular;\n", info, info );
printf( "the solution could not be computed.\n" );
exit( 1 );
}
/* Print solution */
print_matrix( "Solution", n, nrhs, b, ldb );
/* Print details of LU factorization */
print_matrix( "Details of LU factorization", n, n, a, lda );
/* Print pivot indices */
print_int_vector( "Pivot indices", n, ipiv );
exit( 0 );
} /* End of LAPACKE_dgesv Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda ) {
MKL_INT i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
printf( "\n" );
}
}

/* Auxiliary routine: printing a vector of integers */
void print_int_vector( char* desc, MKL_INT n, MKL_INT* a ) {
MKL_INT j;
printf( "\n %s\n", desc );
for( j = 0; j < n; j++ ) printf( " %6i", a );
printf( "\n" );
}
```

another which i created for gnu blas lapack :

```#include<stdio.h>
#include<iostream>
#include "lapacke.h"

using namespace std;

int main()
{
// note, to understand this part take a look in the MAN pages, at section of parameters.
char    TRANS = 'T';
int     INFO=3;
int     LDA = 3;
int     LDB = 3;
int     N = 3;
int     NRHS = 1;
int     IPIV[3] ;

/*	double A[9]=
{
1,2,-1,
2,1,1,
-1,2,1,
};
double B[3]=
{
4,
-2,
2
};
*/

double  A[9] =
{
1,1,1,
1,1,3,
2,1,1
};

double B[3] =
{
1,
2,
3
};

// end of declarations

cout << "compute the LU factorization..." << endl << endl;
//void LAPACK_dgetrf( lapack_int* m, lapack_int* n, double* a, lapack_int* lda, lapack_int* ipiv, lapack_int *info );
LAPACK_dgetrf(&N,&N,A,&LDA,IPIV,&INFO);

// checks INFO, if INFO != 0 something goes wrong, for more information see the MAN page of dgetrf.
if(INFO)
{
cout << "an error occured : "<< INFO << endl << endl;
}else{
cout << "solving the system..."<< endl << endl;
// void LAPACK_dgetrs( char* trans, lapack_int* n, lapack_int* nrhs, const double* a, lapack_int* lda, const lapack_int* ipiv,double* b, lapack_int* ldb, lapack_int *info );
dgetrs_(&TRANS,&N,&NRHS,A,&LDA,IPIV,B,&LDB,&INFO);
printf("IPIV= %d %d %d \n",IPIV[0],IPIV[1],IPIV[2]);
if(INFO)
{
// checks INFO, if INFO != 0 something goes wrong, for more information see the MAN page of dgetrs.
cout << "an error occured : "<< INFO << endl << endl;
}else{
cout << "print the result : {";
int i;
for (i=0;i<N;i++)
{
cout << B << " ";
}
cout << "}" << endl << endl;
}
}

cout << "program terminated." << endl << endl;
return 0;
}```
```compute the LU factorization...
solving the system...
IPIV= 1 3 3
print the result : {2 -1.5 0.5 }
program terminated.```
```outputs are:

LAPACKE_dgesv (row-major, high-level) Example Program Results
Solution
2.00
-1.50
0.50

Details of LU factorization
2.00   1.00   1.00
0.50   0.50   2.50
0.50   1.00  -2.00

Pivot indices
3      2      3

```

can intel's ipiv differ from gnu's ipiv (different algorithm) or, there is some error in my code ?

Regards
Puneet

Black Belt
166 Views

I regret to say that to me the report that you presented was a confusing jumble of source code and results, and I no longer know what it is that you wish to compare for consistency. Taking your thread title as the primary topic, let us take the following short code as our test case.

```#include <stdio.h>

#define N 5
double A[N*N] = {
6.80, -2.11,  5.66,  5.97,  8.23,
-6.05, -3.30,  5.36, -4.44,  1.08,
-0.45,  2.58, -2.70,  0.27,  9.04,
8.32,  2.71,  4.35, -7.17,  2.14,
-9.67, -5.14, -7.26,  6.08, -6.87
};
double b = {4.02,  6.19, -8.22, -7.57, -3.03};

main(){
int piv,info,n,nrhs,i;
n=N; nrhs=1;
dgesv_(&n,&nrhs,A,&n,piv,b,&n,&info);
if(info){
fprintf(stderr,"DGESV returned info=%d\n",info);
exit(1);
}
printf("Solution from DGESV:\n");
for(i=0; i<n; i++)printf("%12.4e ",b);
printf("\nPivots: %3d %3d %3d %3d %3d\n",piv[0],piv[1],piv[2],piv[3],piv[4]);
}
```

After building and running this program using (i) Intel C-32-bit on Windows 10, and (ii) GCC 4.9 and Lapack from Cygwin-64, I obtain identical results (aside from Intel C printing three and GCC printing two digits in the exponent fields).

```s:\lang\mkl>icl /Qmkl xdgesv.c
Intel(R) C++ Intel(R) 64 Compiler for applications running on IA-32, Version 16.0.0.110 Build 20150815

xdgesv.c
Microsoft (R) Incremental Linker Version 12.00.31101.0

-out:xdgesv.exe
-libpath:C:\LANG\INTEL16\compilers_and_libraries_2016.0.110\windows\mkl\lib\ia32_win
xdgesv.obj

s:\lang\mkl>xdgesv
Solution from DGESV:
-8.0071e-001 -6.9524e-001  5.9391e-001  1.3217e+000  5.6576e-001
Pivots:   5   5   3   4   5
---------------
\$ gcc xdgesv.c -llapack -lblas

\$ ./a
Solution from DGESV:
-8.0071e-01  -6.9524e-01   5.9391e-01   1.3217e+00   5.6576e-01
Pivots:   5   5   3   4   5

```