lapack sgetrs tutorial

I have a problem to solving system linear equation Ax=b using lapack. I am planned to use sgetrs subroutine in lapack to solving system linear equation Ax=b.

sgetrs descriptions is a subroutine to Solves a system of linear equations with an LU-factored square coefficient matrix, with multiple right-hand sides. So, we need to process our data with LU-factored before use sgetrs lapack.

I have create a C/C++ code to give tutorial how to call sgetrs lapack from C/C++ code. I have compared the result from this sgetrs lapack with Matlab and give the same result. Please check this code below (call sgetrs lapack from C/C++) :

/*
 * main.c
 *
 *  Created on: 21 Apr, 2016
 *      Author: toto
 */
#include <stdlib.h>
#include <stdio.h>

extern void sgetrf_(int *N, int *NRHS, float *A,
                      int *LDA, int *IPIV, int *INFO );
extern void sgetrs_(char *TRANS, int *N, int *NRHS, float *A,
                      int *LDA, int *IPIV, float *B, int *LDB, int *INFO );

void print1float(char *ca, float *a, int n)
{
	int i;
	fprintf(stderr, "%s:\n", ca);
	for(i=0; i<n; i++)
		fprintf(stderr, "%f \n", a[i]);
	fprintf(stderr, "\n");
}

void print2float(char *ca, float *a, int nrow, int ncol)
{
	int i, j, ii;

	ii = 0;
	fprintf(stderr, "%s:\n", ca);
	for(i=0; i<nrow; i++)
	{
		for(j=0; j<ncol; j++)
		{
			fprintf(stderr, "%f \t", a[ii]);
			ii = ii+1;
		}
		fprintf(stderr, "\n");
	}
	fprintf(stderr, "\n");
}

int main(int argc, char **argv)
{
	float A[16] = {
			0.5377,    0.3188,    3.5784,    0.7254,
		    1.8339,   -1.3077,    2.7694,   -0.0631,
		    -2.2588,   -0.4336,   -1.3499,    0.7147,
		    0.8622,    0.3426,    3.0349,   -0.2050};
	float b[4] = {1,2,3,4};

	int i;
	char trans = 'T';
	int dim = 4;
	int nrhs = 1;
	int LDA = dim;
	int LDB = dim;
	int info;
	int ipiv[4];

	//show input array
	print2float("Input A", A, dim, LDA);
	print1float("Input b" , b, dim);

	//compute LU decomposition
	sgetrf_(&dim, &dim, A, &LDA, ipiv, &info);
	print2float("LU A", A, dim, LDA);

	//solve Ax = b
	sgetrs_(&trans, &dim, &nrhs, A, &LDA, ipiv, b, &LDB, &info);
	print1float("x" , b, dim);

	return(0);
}

We can compile this code (call sgetrs lapack from C/C++) using command :

gcc main.c -lm -llapack -lblas -o main

Sample output after running this code (call sgetrs lapack from C/C++) :

Input A:
0.537700 	0.318800 	3.578400 	0.725400 	
1.833900 	-1.307700 	2.769400 	-0.063100 	
-2.258800 	-0.433600 	-1.349900 	0.714700 	
0.862200 	0.342600 	3.034900 	-0.205000 	

Input b:
1.000000 
2.000000 
3.000000 
4.000000 

LU A:
3.578400 	0.089090 	0.150263 	0.202716 	
2.769400 	-1.554426 	-0.912081 	0.401758 	
-1.349900 	-0.313337 	-2.341749 	-0.475812 	
3.034900 	0.072220 	0.472039 	-0.624637 	

x:
-4.261576 
-2.205967 
2.362486 
-6.147243