Solving Matrices in C using Lapack

Question

So I have two matrices; One is 2 by 2 and the other is 2 by 1. I want to use the external library lapack to solve the linear system and it seems like I need to call the function dgesv_() which is

   Ax = B

And I have to solve for x.

So I am really confused on the notation they are saying to call the function and how to translate it into the two arrays I have now.

 #include <stdlib.h>
 #include <stdio.h>


static double Angstroms[2];
static double Energy[2];
static double ax[2][2];

void file_input ();
void polynomial ();


int main () {

     file_input ();
     polynomial ();
     return 0;
}

 void file_input () {

     float a, b;
     int i;

     FILE * in_file = fopen("H2Mini.txt", "r");
     if (outfile == NULL) {
         printf ("Error file does not exist");
         exit (-1);
     }
     for (i = 0; i <= 1; i++) {
          fscanf(in_file, "%f %f\n", &a, &b);
          Angstroms[i] = a;
          Energy [i] = b;
     }
     fclose(in_file);

}

 void polynomial () {

      int i;
      FILE * outfile = fopen("PolyTest2.txt", "w");
      if (outfile == NULL) {
         printf ("Error file does not exist");
         exit (-1);
      }
     for (i = 0; i <= 1; i++) {
         ax[i][0] = 1;
         fprintf (outfile, "%.8f ", ax[i][0]);

     }
     fprintf (outfile, "\n");
     for (i = 0; i <= 1; i ++) {
         ax[i][1] = Angstroms[i];
         fprintf (outfile, "%.8f ", ax[i][1]);
     }
  }

So I have my in File is this

    [2.00000000 3.00000000
     6.00000000 5.00000000]

The ax array is going to look like this

    [1.00000000 1.00000000
     2.00000000 6.00000000]

And the Energy array is this

    [3.00000000 5.00000000]

My ax array is the Ax term in the Ax = b equation and my b term is the Energy array.

I did look over their function documentation and it is a little confusing in implementing it. In dgesv (n, nrhs, a, lda, ipiv, b, ldb, info)

Any really clear examples of this code would be super helpful!

"So I have made two matrices; One is **2 by 2** and the other is **1 by 1**." - and you want to multiply these? — user2357112, May 05 '15 at 05:33
You find the C LAPACK interface (from Netlib) easier to use. — Jeff Hammond, May 06 '15 at 19:47
Maybe... The CLAPACK library was built using a Fortran to C conversion utility. It is not without pitfalls. You are far better served managing your own matrix multiplications in c without CLAPACK, unless you get to the point your needs are so complex a pre-written library, such as CLAPACK becomes an option. — David C. Rankin, May 11 '15 at 16:29

score 0 · Answer 1 · answered May 11 '15 at 14:37

I am not sure what you mean by 'The ax array is going to look like this'? Ax should be a vector shouldn't it?

However in general the call to dgesv using the C wrapper LAPACKE should look like this:

LAPACKE_dgesv(
            LAPACK_ROW_MAJOR, // The storage format you use, usually row major in c.
            n,                // Dimensions of A (nxn)
            1,                // Number of b's, here always one because you want to solve for a single right hand side.
            A,                // The input matrix A, it will be destroyed in the call so best make a copy.
            n,                // LDA basically the length of each column of A in memory. Normally this is also n if A is nxn.
            p,                // An additional array where lapacke saves pivot ordering. Probably not of any interest to you. Just give it an array of size n. 
            bToX,             // On input input the rhs b, on output this will be the solution vector x
            1 );

Solving Matrices in C using Lapack

1 Answers1