- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
i have two codes for solving linear equations, one from intel:
/*******************************************************************************
* Copyright (C) 2009-2015 Intel Corporation. All Rights Reserved.
* 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
* permission. No license under any patent, copyright or other intellectual
* 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 ?
Awaiting your reply
Regards
Puneet
Link Copied
1 Reply
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 Copyright (C) 1985-2015 Intel Corporation. All rights reserved. xdgesv.c Microsoft (R) Incremental Linker Version 12.00.31101.0 Copyright (C) Microsoft Corporation. All rights reserved. -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
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page