LAPACK C (mkl) dptsv row major/column major: does it make a difference for vectors?

Question

Folks, I'm calling LAPACKE_dptsv

If all arguments are 1d arrays, does it matter if I call the data LAPACK_ROW_MAJOR or LAPACK_COL_MAJOR ?

Answer 1

The function LAPACKE_dptsv() corresponds to the lapack function dptsv(), which does not feature the switch between LAPACK_ROW_MAJOR and LAPACK_COL_MAJOR. dptsv() is implemented for column-major ordering, corresponding to matrices in Fortran, while most of C matrices are row-major. So LAPACKE_dptsv(LAPACK_ROW_MAJOR,...) performs the following steps :

• transpose the right-end side b
• call dptsv() of Lapack
• transpose the response (that is b again)

You can check this in the source of Lapacke, in /lapacke/src/lapacke_dptsv_work.c.

There is one question left : *does it largely affect the wall-clock time ? * Looking at dpttrs_8f_source, it is possible : the decomposition L*D*L**T is performed by a single for loop (+loop unrolling). A piece of code is therefore necessary to answer the question. The following code is compiled by gcc main.c -o main -llapacke -llapack -lblas

#include <stdio.h>
#include "lapacke.h"
#include <malloc.h>
#include <time.h>

int main ()
{
    //double a[3][2] = {{1,0},{1,1},{1,2}};

    double **outputArray;
    int designs=3;
    int i,j;
    lapack_int info,n,ldb,nrhs;

    n = 420000;
    nrhs = 42;

    //double outputArray[3][1] = {{6},{0},{0}};
    double* ad=malloc(n*sizeof(double));
    if(ad==NULL){printf("malloc failed\n");exit(1);}
    for(i=0;i<n;i++){
        ad[i]=3;
    }
    double* ae=malloc((n-1)*sizeof(double));
    if(ae==NULL){printf("malloc failed\n");exit(1);}
    for(i=0;i<n-1;i++){
        ae[i]=-1;
    }

    double* b=malloc(n*nrhs* sizeof(double));
    if(b==NULL){printf("malloc failed\n");exit(1);}

    for(j=0;j<nrhs;j++){
        for(i=0;i<n;i++){
            b[i*nrhs+j]=i+2*j;
        }
    }

    ldb=nrhs;

    clock_t t;
    t = clock();
    info = LAPACKE_dptsv(LAPACK_ROW_MAJOR,n,nrhs,ad,ae,b,ldb);
    if(info!=0){printf("failed, info %d\n",info);}
    t = clock() - t;
    printf ("LAPACK_ROW_MAJOR :  %d clicks (%f seconds).\n",t,((float)t)/CLOCKS_PER_SEC);

    for(i=0;i<n;i++){
        ad[i]=3;
    }
    if(ae==NULL){printf("malloc failed\n");exit(1);}
    for(i=0;i<n-1;i++){
        ae[i]=-1;
    }

    double* b2=malloc(n*nrhs* sizeof(double));
    if(b2==NULL){printf("malloc failed\n");exit(1);}
    for(j=0;j<nrhs;j++){
        for(i=0;i<n;i++){
            b2[j*n+i]=i+2*j;
        }
    }

    t = clock();
    ldb=n;
    info = LAPACKE_dptsv(LAPACK_COL_MAJOR,n,nrhs,ad,ae,b2,ldb);
    if(info!=0){printf("failed, info %d\n",info);}
    t = clock() - t;
    printf ("LAPACK_COL_MAJOR :  %d clicks (%f seconds).\n",t,((float)t)/CLOCKS_PER_SEC);

    double delta=0,temp,deltal=0;
    for(i=0;i<n;i++){
        deltal=0;
        for(j=0;j<nrhs;j++){
            temp=(b[i*nrhs+j]-b2[j*n+i]);
            deltal+=temp*temp;
        }
        delta+=deltal;
    }
    printf("delta %g\n",delta);
    free(ad);
    free(ae);
    free(b);
    free(b2);
    return (info);
} 

My output is :

LAPACK_ROW_MAJOR : 770000 clicks (0.770000 seconds).

LAPACK_COL_MAJOR : 310000 clicks (0.310000 seconds).

So LAPACKE_dptsv() runs almost twice faster with LAPACK_COL_MAJOR with `nbrhs=42

If the number of rhs is reduced to one (and n larger) :

LAPACK_ROW_MAJOR : 250000 clicks (0.250000 seconds).

LAPACK_COL_MAJOR : 180000 clicks (0.180000 seconds).

LAPACK_COL_MAJOR and LAPACK_ROW_MAJOR lead to approximately the same wall-clock time with a single RHS. And the output is the same.

I have not installed intel mkl on my pc and i am curious about the way it changes the conclusion of this test...