Using DGESV in Lapacke [solved]

[ctcsback]

As requested by Julie, I'm reposting my previous problem using DGESV in lapacke.

For some reason, my previous implementation wasn't giving understandable results in lapacke, after having been running it correctly under lapack. Using Matlab as a check, I'd try running row major implementation of dgesv and get garbage, even though column major settings would check off with Matlab. Needless to say, it was frustrating. Attached is my older code

Code: Select all

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
//#include <cblas.h>
#include <lapacke.h>

double A[36] =
{
4, 3, 0, 4, 0, 3,
2, 3, 3, 2, 2, 0,
0, 3, 2, 2, 0, 3,
3, 3, 2, 3, 0, 0,
1, 2, 3, 4, 0, 2,
2, 4, 2, 4, 3, 2
};

double C[24] =
{
34, 36, 33, 39,
33, 18, 26, 23,
28, 24, 17, 23,
30, 22, 29, 21,
31, 32, 24, 28,
44, 32, 35, 44
};

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

void eq_solve(double *a, double *b, int n, int m)
{
  int info, i;
  int lda, ldb, *ipiv;
  lda = n;
  ldb = n;
  ipiv = malloc(sizeof(int)*n);

  info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, a, lda, ipiv, b, ldb);

  if (info != 0)
  {
    fprintf(stderr, "dgesv: info = %d\n", info);
  }
  assert(info==0);

  free(ipiv);
}

int main(void)
{
  eq_solve(A,C,6,4);

// Print and compare B and C here

  return 0;
}


This was compiled with the lapacke.h header and the libraries in folder under command: gcc test.c -I../include lapacke.a -llblas -llapack

The garbage answer that I kept getting was:

14.2353, 14.1176, 9.17647, 14.4118,
33, 18, 5.48416, 4.97285,
14.8824, 5.90498, 17, 23,
8.85068, 6.88688, 4.17647, 6.10407,
31, 32, -17.6199, -14.3484,
-12.1765, -13.7195, 35, 44,

After extensive research, I stumbled onto a file from intel's MKL website.
http://software.intel.com/sites/product ... _row.c.htm

I was doing a direct translation of my old Lapack call to the Lapacke call without considering col-major to row-major. LDB is no longer equal to n. Instead it is nrhs now. *facepalm*

Aside from that, I spent a good deal of timing trying to get both Lapacke and Cblas working for my C code so that all my matrices can now be supported in row major form. No more pre and post processing of transposing matrices :).

I got it to work, but not without a bunch of initial errors. For some reason, my makefile MUST link the libraries in a specific order - that is :

Code: Select all

LDLIBS := -lm -lgfortran -lmetis -lpthread -lcblas_LINUX -llapacke_LINUX -llapack_LINUX -lblas_LINUX

Switching the order always results in some undefined reference to dgemm_ or dgetrf_ or something like this: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1669 .

I hope this helps.

[admin]

I was doing a direct translation of my old Lapack call to the Lapacke call without considering col-major to row-major. LDB is no longer equal to n. Instead it is nrhs now. *facepalm*

Exactly that was the problem.
Here is a quick example code that works in Colomn Major first then in Row Major - Both are equivalent and give the same output.

Code: Select all

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <assert.h>
#include <time.h>
#include "lapacke.h"

double A[36] =
{
4, 3, 0, 4, 0, 3,
2, 3, 3, 2, 2, 0,
0, 3, 2, 2, 0, 3,
3, 3, 2, 3, 0, 0,
1, 2, 3, 4, 0, 2,
2, 4, 2, 4, 3, 2
};

double Arow[36] =
{
     4,     2,     0,     3,     1,     2,
     3,     3,     3,     3,     2,     4,
     0,     3,     2,     2,     3,     2,
     4,     2,     2,     3,     4,     4,
     0,     2,     0,     0,     0,     3,
     3,     0,     3,     0,     2,     2,
};


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

double Brow[24] =
{
     2,     3,     2,     1,
     2,     1,     4,     4,
     3,     3,     3,     2,
     2,     2,     4,     4,
     4,     1,     2,     0,
     0 ,    0,     0,     4,
};

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

        else {
        printf("m=%d,n=%d\n",m,n);
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", a[j+i*lda] );
                printf( "\n" );
        }
        }
       
}

void eq_solve(double *a, double *b, int n, int nrhs, char major)
{
  int info, i;
  //char trans = 'n';
  lapack_int lda, ldb, *ipiv;
  lda = n;
  ipiv = malloc(sizeof(int)*n);
        /* Print A */
  print_matrix( "A", major, n, n, a, lda );

  if (major == 'C'){
  ldb = n;
        /* Print B */
  print_matrix( "B", major, n, nrhs, b, ldb );
  //EQUIVALENT STATEMENT
  info = LAPACKE_dgesv(LAPACK_COL_MAJOR, (lapack_int) n, (lapack_int) nrhs, a, lda, ipiv, b, ldb);
  //dgesv_( &n, &nrhs, a, &lda, ipiv, b, &ldb, &info );
  }
  else {
  ldb = nrhs;
        /* Print B */
  print_matrix( "B", major, n, nrhs, b, ldb );
  info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, (lapack_int) n, (lapack_int) nrhs, a, lda, ipiv, b, ldb);
   }
 

  if (info != 0)
  {
    fprintf(stderr, "dgesv: info = %d\n", info);
  }
  assert(info==0);
  printf("dgesv passed\n");
  //debug_print("b using lapacke", b,n,nrhs,0);
  //print answer here
  print_matrix( "Answer", major, n, nrhs, b, ldb );

  free(ipiv);
}

int main(void)
{
  printf("==============================\n");
  printf("\tSOLVING IN COL MAJOR\n");
  printf("==============================\n");
  eq_solve(A,B,6,4,'C');
  printf("\n==============================\n");
  printf("\tSOLVING IN ROW MAJOR\n");
  printf("==============================\n");
  eq_solve(Arow,Brow,6,4,'R');

  return 0;
}




To compile it

Code: Select all

gcc test.c -Imydirlapacke/include  -Lmydirlib -lliblapacke -lliblapack -lllibblas -lgfortran


In the mydirlapacke/include, I have the lapack.h
In mydirlib, I have liblapacke.a , libblas.a and liblapack.a

[Julien Langou]

Hi,

Thanks for the feedback! I am sure the post will show to be useful by others!

1) When switching from col major to row major, the leading dimension changes. Be aware!

2) For the ordering of the library, on our side, the constrains are:
you need lapacke in front of lapack
you need lapack in front of blas
you need cblas in front of blas
(as you pointed out)

Julien.

[ctcsback]

I don't know exactly why, but my program has slowed down significantly. Before using Lapacke, I'd have to solve my matrix problem and then do some post-processing to transpose my matrix back to row major, which I believe to have taken a large chunk of my time. Is Lapacke supposed to be slower - or could it be due to other factors?

For instance, I'm now linking all my lapack libraries in a local directory, whereas before i was just using -llapack -lblas which was getting the libraries from /usr/lib/.

For a given problem, when running on a single thread, my previous implementation is faster by a factor of two, and it gets worse as the number of threads increase (ran on 4 and was almost 3x faster).

[admin]

There is no reason to notice a drop in performance when switching to LAPACKE.
I believe that now you are using the Reference BLAS and this is what is slowing you down. The BLAS from /usr/lib should be an optimized one.
You should try to link with your LAPACK and BLAS from /usr/lib/.

Please let me know how it goes.
This is very interesting to us.
Julie

[ctcsback]

Sorry, I just realize my runs were on a different machines. It is interesting, since the faster runs are on a older and slower AMD server architecture while the new runs were on my new laptop running Intel's latest architecture. Running my old code on my laptop now, the numbers seem more reasonable ~ 25% speedup without the extra transposing using lapacke.

I think you are right about the /usr/lib/ libraries being optimized. For my given program, when I change -lblas_LINUX to -lblas, the runtime decreases from 25.81 to 12.93, and then changing -llapack_LINUX to -llapack, decreased runtime even more to 11.27. Should ATLAS give even faster results then?

[Julien Langou]

Hi,

1) It is highly probable that the optimized BLAS that is shipped with your Linux distribution is coming from ATLAS.
This is kind of the only open source optimized BLAS library around.
(OK, there is gotoblas, but I think gotoblas license does not allow a linux distribution to ship it.)

2) To be sure try something like "nm libblas.a | grep -i atlas"

3) if libblas.a is an atlas library then ... some would argue that you still want to do your own atlas install.
Atlas is an autotuning system, so the blas you get is optimized for your very own specific computer.
Well, I would personally would go with whatever is shipped within the Linux distribution and do not do my own install.
You can try if you want.

4) If it's not atlas, well, I'd like to know what it is.

5) Yep the lapack_LINUX is the one from netlib. When do we remove this _LINUX from the name of the library .... ????

Julien