Noob: Lapack API - using methods

Open discussion regarding features, bugs, issues, vendors, etc.

Noob: Lapack API - using methods

Postby Mat » Sat Aug 19, 2006 10:10 am

Hello,
i'm new to Lapack and to this forum.

I'm sorry for bothering with an easy question i'm not able to answer
myself. I'm searching for something like the API for Lapack. For example i would like to compute a simple QR factorization but i don't know how to implement it in my C++ code. I found the User's Guide here but it was not helpful at all.
It would be very kind if somebody could help me. Perhaps i'm just to blind to see some little detail.

Thanks a lot
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julie » Sat Aug 19, 2006 12:25 pm

Hi Mat,

please take a look at the LAPACK FAQ http://www.netlib.org/lapack/faq.html#1.1 to see you you don't have the LAPACK library already installed on you machine.

The routine to perform a QR factorization of a general matrix is DGEQRF ( D for Double precision solve).
Here is the documentation for the Fortran routine: http://www.netlib.org/lapack/double/dgeqrf.f
Here is an example from the "IBM Lapack": ESSL http://www.cepba.upc.edu/docs/ibm_doc/manessl/essl368.html

If you are looking for a C++ Wrapper for Lapack routine, you can have a look at http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=147&sid=d634dcd28d862f9d1066bfa84b79da9a

But you can directly call Fortran routines from C++. Here is how:

The extern "C" directive

First of all, the extern "C" linkage specification must be used to declare any modules that are not written in C++.
So you will have to define the Fortran routine in your C++code:
Code: Select all
extern "C" void dgeqrf_ (const int &m, const int &n, double da[], const int &lda, double dtau[], double dwork[], const int &ldwork, int &info)

Case sensitivity
C++ is a case sensitive language, Fortran is not. The name of the Fortran routine is mangled by the compiler. All references to Fortran symbols may be specified in lowercase or in uppercase (depends on the compiler) in C++ calls.

The underscore (_)

In many Fortran compilers, an underscore (_) is appended by default at the end of definitions and references to externally visible symbols (subroutines, functions, etc.). This happens for Solaris and Digital UNIX f77 compilers, for g77 compiler and for hepf77 compiler interface; on HP-UX the +ppu option must be added to the f77 compiler to activate this feature; on AIX -qextname must be added to the xlf compiler. The appended underscore is a standard for HEP Fortran libraries (e.g. CERNLIB). In a C++ call the Fortran subroutines or functions must be specified with an appended underscore.

Passing parameters
An important thing to remember is that C++ passes all parameters by value (except arrays and structures). Fortran passes them by reference. It is therefore necessary to specify in the C++ function prototypes that the Fortran subroutines and functions expect call-by-reference arguments using the address-of (&) operator.

Arrays
C++ stores arrays in row-major order, whereas Fortran stores arrays in column-major order. The lower bound for C++ is 0, for Fortran is 1. This mean that the Fortran array element myarray(1,1) corresponds to the C++ array element myarray[0][0]; whereas myarray(6,8) corresponds to myarray[7][5];

How to compile and link

Obviously the C++ and Fortran source files must be compiled separately using the appropriate compiler. Linker must be called via the C++ compiler as the main program lies in the C++ source files:
Code: Select all
g77 -c myfortran.f
g++ -o myprogram myprogram.cc myfortran.o -lg2c  -lm

and for a library:
Code: Select all
g++ -o myprogram myprogram.cc PATH_TO_LIBRARY/lapack.a -lg2c -lm
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

Postby Mat » Sun Aug 20, 2006 2:28 pm

:D
Thanks very much for the detailed description. I'm really happy somebody answered in so detailed help!

I'm looking forward to check it out!
Thanks
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Mat » Sun Aug 20, 2006 3:42 pm

Its me again :cry:

sorry but i am in need of some further help. I really read the documentation but because of my english (which isn't very good) i have
still problems.
Now i know how to implement the Lapack routines but i have great problems in compiling the code correctly.

My programm looks this like:
Code: Select all
extern "C" void dgeqrf_(const int &m, const int &n, double da[], const int &lda, double dtau[], double dwork[], const int &ldwork, int &info);

Class::Class()
{
int n = 2,
       m = 6,
       lda = 6,
       info;
   double work[5],
         tau[5],
          A[25];
       
       A[0] = 0.0;
       A[1] = 2.0;
       A[2] = 2.0;
       A[3] = -1.0;
       A[4] = 2.0;
       A[5] = -1.0;
       A[6] = 0.0;
       A[7] = 1.5;
       A[8] = 2.0;
       A[9] = -1.0;
       A[10] = 2.0;
       A[11] = -1.0;
       

           
   std::cout << "input A:" << std::endl;
   for(int i=0; i < m; i++){
      for(int j=0; j < n; j++){
          std::cout << A[i*n +j] << std::endl;
      }
      std::cout << ""<< std::endl;
   }
   std::cout << "------------------------------------------------------------\n"<< std::endl;

   
   dgeqrf_(m,n,A,lda,tau,work,n,info);
   
   
   
   std::cout << "output A:"<< std::endl;
   for(int i=0; i < m; i++){
       for(int j=0; j < n; j++){
           std::cout << A[i*n+j] << std::endl;
       }
       std::cout << "" << std::endl;
   }
   std::cout << "------------------------------------------------------------\n"<< std::endl;
}


And i use this compiler option:

Code: Select all
g++ -c -O3 -lg2c -lm -o myProg.o myProg.cpp


The error message is:
Code: Select all
g++: -lg2c: input files of the binder are unused, because there is no binding
g++: -lm:  input files of the binder are unused, because there is no binding

myProg.cpp: undefined reference to `dgeqrf_'


In the description above in section: "How to compile and link"
there is :
Code: Select all
g77 -c myfortran.f

Can you explain me which fortran file this should be? Do i have to write my own fortran file?

Addiotionally i don't know what you mean with library
Code: Select all
g++ -o myprogram myprogram.cc PATH_TO_LIBRARY/lapack.a -lg2c -lm

and in my Lapack folder there is no lapack.a file :(

And do i need the C++ wrapper? I'm sorry but i'm not familiar with wrappers - what is this and for what do i need it?

Thanks for patient help :)
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julie » Sun Aug 20, 2006 4:33 pm

You are not far....
Hers is the declaration...
Code: Select all
extern "C" void dgeqrf_(int *m, int *n, double da[], int *lda, double dtau[], double dwork[], int *ldwork, int *info);

here is the call
Code: Select all
dgeqrf_(&m,&n,A,&lda,tau,work,&n,&info);

here is the compilation command
Code: Select all
g++ -O3 -lg2c -lm -o exe myProg.cpp ~/lib/liblapack.a ~/lib/libblas.a

You NEED the LAPACK and BLAS librairies...what do you have inside you LAPACK folder?
(~/lib is the path on my machine to my LAPACK and BLAS libraries)
Julie
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

Postby Mat » Sun Aug 20, 2006 4:58 pm

Thanks again julie,

i found my liblapack.a and libblas.a files. They are under /usr/lib/atlas.

The compilation command looks now:
Code: Select all
g++ -O3 -lg2c -lm -o result myProg.cpp /usr/lib/atlas/liblapack.a /usr/lib/atlas/libblas.a


But i receive the following error which i don't know to handle with:
Code: Select all
/usr/lib/gcc-lib/i486-linux/3.3.5/../../../crt1.o(.text+0x18): In function `_start':
../sysdeps/i386/elf/start.S:98: undefined reference to `main'
collect2: ld returned 1 exit status


It looks to become problematic :(
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julie » Sun Aug 20, 2006 5:03 pm

Mat,

exe is the name I gave to my executable program.
Feel free to call your program whatever you like.
Here is the compilation command you should use:

Code: Select all
g++ -O3 -lg2c -lm -o WhatEverProgramName myProg.cpp -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas


You have to use the LAPACK library from DEBIAN (in /usr/lib), not the one from Atlas (they only provide part of the library).
For information, your LAPACK library is coming from the debian RPM, thanks to Camm Maguire.

Here is my complete program.
I manage to compile and run it.
Code: Select all
#include <iostream>
extern "C" void dgeqrf_(int *m, int *n, double da[], int *lda, double dtau[], double dwork[], int *ldwork, int *info);
main{
   int n = 2,
       m = 6,
       lda = 6,
       info;
   double work[5],
         tau[5],
          A[25];

       A[0] = 0.0;
       A[1] = 2.0;
       A[2] = 2.0;
       A[3] = -1.0;
       A[4] = 2.0;
       A[5] = -1.0;
       A[6] = 0.0;
       A[7] = 1.5;
       A[8] = 2.0;
       A[9] = -1.0;
       A[10] = 2.0;
       A[11] = -1.0;



   std::cout << "input A:" << std::endl;
   for(int i=0; i < m; i++){
      for(int j=0; j < n; j++){
          std::cout << A[i*n +j] << std::endl;
      }
      std::cout << ""<< std::endl;
   }
   std::cout << "------------------------------------------------------------\n"<< std::endl;


   dgeqrf_(&m,&n,A,&lda,tau,work,&n,&info);


   if (info == 0)
   {
     std::cout << "output A:"<< std::endl;
     for(int i=0; i < m; i++){
       for(int j=0; j < n; j++){
           std::cout << A[i*n+j] << std::endl;
       }
       std::cout << "" << std::endl;
     }
   }
   else
           std::cout << "Problem!!! info = " << info << "\n";
   std::cout << "------------------------------------------------------------\n"<< std::endl;
}



Julie
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

Postby Mat » Sun Aug 20, 2006 5:16 pm

Thanks a lot. Yes i'm very confused at the moment sorry.

The example you wrote works fine:
Code: Select all
#include <iostream>
extern "C" void dgeqrf_(int *m, int *n, double da[], int *lda, double dtau[], double dwork[], int *ldwork, int *info);
main{
   int n = 2,
       m = 6,
       lda = 6,
       info;
   double work[5],
         tau[5],
          A[25];
 
       A[0] = 0.0;
       A[1] = 2.0;
       A[2] = 2.0;
       A[3] = -1.0;
       A[4] = 2.0;
       A[5] = -1.0;
       A[6] = 0.0;
       A[7] = 1.5;
       A[8] = 2.0;
       A[9] = -1.0;
       A[10] = 2.0;
       A[11] = -1.0;
 
 
 
   std::cout << "input A:" << std::endl;
   for(int i=0; i < m; i++){
      for(int j=0; j < n; j++){
          std::cout << A[i*n +j] << std::endl;
      }
      std::cout << ""<< std::endl;
   }
   std::cout << "------------------------------------------------------------\n"<< std::endl;
 
 
   dgeqrf_(&m,&n,A,&lda,tau,work,&n,&info);
 
 
   if (info == 0)
   {
     std::cout << "output A:"<< std::endl;
     for(int i=0; i < m; i++){
       for(int j=0; j < n; j++){
           std::cout << A[i*n+j] << std::endl;
       }
       std::cout << "" << std::endl;
     }
   }
   else
           std::cout << "Problem!!! info = " << info << "\n";
   std::cout << "------------------------------------------------------------\n"<< std::endl;
}



BUT the following example does not work:
Code: Select all
#include <iostream>
#include "myProg.h"
 extern "C" void dgeqrf_(int *m, int *n, double da[], int *lda, double dtau[], double dwork[], int *ldwork, int *info);
MyProg::MyProg()
{
   
    int n = 2,
        m = 6,
        lda = 6,
        info;
    double work[5],
          tau[5],
           A[25];
 
        A[0] = 0.0;
        A[1] = 2.0;
        A[2] = 2.0;
        A[3] = -1.0;
        A[4] = 2.0;
        A[5] = -1.0;
        A[6] = 0.0;
        A[7] = 1.5;
        A[8] = 2.0;
        A[9] = -1.0;
        A[10] = 2.0;
        A[11] = -1.0;
 
 
 
    std::cout << "input A:" << std::endl;
    for(int i=0; i < m; i++){
       for(int j=0; j < n; j++){
           std::cout << A[i*n +j] << std::endl;
       }
       std::cout << ""<< std::endl;
    }
    std::cout << "------------------------------------------------------------\n"<< std::endl;
 
 
    dgeqrf_(&m,&n,A,&lda,tau,work,&n,&info);
 
 
    if (info == 0)
    {
      std::cout << "output A:"<< std::endl;
      for(int i=0; i < m; i++){
        for(int j=0; j < n; j++){
            std::cout << A[i*n+j] << std::endl;
        }
        std::cout << "" << std::endl;
      }
    }
    else
            std::cout << "Problem!!! info = " << info << "\n";
    std::cout << "------------------------------------------------------------\n"<< std::endl;
 }
 
MyProg::~MyProg() {}

with the header:

Code: Select all

class MyProg
{
 
    public:
       
        MyProg();
       
        ~MyProg();
           
};


the error message is as already posted:
Code: Select all
/usr/lib/gcc-lib/i486-linux/3.3.5/../../../crt1.o(.text+0x18): In function `_start':
../sysdeps/i386/elf/start.S:98: undefined reference to `main'
collect2: ld returned 1 exit status



Do you know why?
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julie » Sun Aug 20, 2006 8:05 pm

Mat,

a C++ program needs a main part....
what you did is write a class...you need a ?"main" to call your class.

Please "google" for some C++ tutorial, they will tell you everything about C++ programming.
Maybe switch to C, it may be easier. depending on what you want to do, to start programming.

My advice is for you to take some time reading a good programming book.
It will really worth it!
Julie
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

Postby Mat » Mon Aug 21, 2006 1:06 pm

Hello Julie,

thanks for your reply.
I knows but the code i posted was just a little bit of my project i work on at the moment. It would be to much to post it in here. Surely my project has a main method but there i get the error message:
Code: Select all
g++ -c -O3 -lg2c -lm -o UIMain.o UIMain.cpp -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas
g++: -lg2c: linker input file unused because linking not done
g++: -lm: linker input file unused because linking not done
g++: -llapack: linker input file unused because linking not done
g++: -lf77blas: linker input file unused because linking not done
g++: -latlas: linker input file unused because linking not done
g++ -c -O3 -lg2c -lm -o MyProg.o MyProg.cpp -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas
In file included from MyProg.cpp:2:
MyProg.h:17:3: warning: no newline at end of file
MyProg.cpp:66:2: warning: no newline at end of file
g++: -lg2c: linker input file unused because linking not done
g++: -lm: linker input file unused because linking not done
g++: -llapack: linker input file unused because linking not done
g++: -lf77blas: linker input file unused because linking not done
g++: -latlas: linker input file unused because linking not done
MyProg.o(.text+0x1de): In function `MyProg::MyProg[not-in-charge]()':
MyProg.cpp:43: undefined reference to `dgeqrf_'
MyProg.o(.text+0x48e): In function `MyProg::MyProg[in-charge]()':
MyProg.cpp:43: undefined reference to `dgeqrf_'
collect2: ld returned 1 exit status
make[1]: *** Error 1


But the file i implemented in my project which shall use the QR factorization has completely the form of MyProg as already posted.
AND there is an class instantiation in another class. Sorry for the misunderstandings and the wrong post yesterday.
My Makefile looks a bit more confusing so the compilation command is this like:


Code: Select all
.cpp.o:
   $(CXX) -c -O3 -lg2c -lm $(CXXFLAGS) $(INCPATH) -o $@ $< -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas

.cc.o:
   $(CXX) -c -O3 -lg2c -lm $(CXXFLAGS) $(INCPATH) -o $@ $< -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas

.cxx.o:
   $(CXX) -c -O3 -lg2c -lm $(CXXFLAGS) $(INCPATH) -o $@ $< -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas

.C.o:
   $(CXX) -c -O3 -lg2c -lm $(CXXFLAGS) $(INCPATH) -o $@ $< -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas

.c.o:
   $(CC) -c -O3 -lg2c -lm $(CFLAGS) $(INCPATH) -o $@ $< -L/usr/lib/ -llapack -L/usr/lib/atlas -lf77blas -latlas


Whats wrong with it?

P.S. CXX = g++
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Postby Julie » Mon Aug 21, 2006 1:39 pm

Hi Mat

Your problem is that you have to do first the compilations of all your files and then the linking.
Here are my main program, my compilation command, my linking sequence and the result of the execution.
Adapt them to your project.

Main.cpp
Code: Select all
#include "myProg.h"
main()
{
        MyProg julie;
}

Compilations comands
Code: Select all
g++ -c myProg.cpp
g++ -c Main.cpp

Linking sequence
Code: Select all
g++ -o MyProgram  Main.o myProg.o ~/lib/liblapack.a ~/lib/libblas.a -lg2c                                     

Execution
Code: Select all
./MyProgram
input A:
0
2

2
-1

2
-1

0
1.5

2
-1

2
-1

------------------------------------------------------------

output A:
-3.74166
0.534522

0.534522
-0.267261

0.534522
-0.267261

-3.4744
0.422577

-0.183216
0.091608

-0.183216
0.091608

------------------------------------------------------------


Hope it helps
Julie
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

Postby Mat » Mon Aug 21, 2006 2:25 pm

:D
Thanks Julie! I had to adapt my big Makefile at the right place. The LDFLAGS had to be changed so now it works all fine! I get the same result for the QR factorization as you.

I'm very thankful for your help! Perhaps someday i can help you - who knows :P
Mat
 
Posts: 47
Joined: Sat Aug 19, 2006 9:54 am

Why got different refult from IBM doc.

Postby wxiong » Mon Aug 28, 2006 3:37 pm

Hi, julie:

I am wei. Found you and this lapack development forum from the google search.

I am trying the dgeqrf funtion. I have got exactly same output A as yours. However, I am wondering but it is different from the example from the "IBM lapack" even with the same input. http://www.cepba.upc.edu/docs/ibm_doc/m ... sl368.html .

In addition,I kind of know how to get matrix Q but not R in A = Q * R, do you know how to get matrix R in A = Q*R?

Thanks for your help in advance.

Wei.
wxiong
 
Posts: 3
Joined: Mon Aug 28, 2006 3:04 pm

Postby Julie » Mon Aug 28, 2006 4:44 pm

Hello Wei,

the program that was derived previously is not exactly solving the problem given as an example by ESSL ( http://www.cepba.upc.edu/docs/ibm_doc/m ... sl368.html ).

Although we started from ESSL's example, I took Mat code and the goal was to help him to have his compilation right.

If you want to initialize the matrix to reproduce ESSL example, you will need to write:
Code: Select all
// first column
       A[0] =   0.0;
       A[1] =   2.0;
       A[2] =   2.0;
       A[3] =   0.0;
       A[4] =   2.0;
       A[5] =   2.0;
// second column
       A[6] =   2.0;
       A[7] =  -1.0;
       A[8] =  -1.0;
       A[9] =   1.5;
       A[10] = -1.0;
       A[11] = -1.0;

Notice that you input A by column since A is to be given to DGEQRF a Fortran routine (so column-major format).

Then when you print the coefficient of the matrix A (before and after the call to DGEQRF), you might want to write:
Code: Select all
     for(int i=0; i < m; i++){
       for(int j=0; j < n; j++){
         std::cout << A[i +j*m] << "\t" ;
       }
       std::cout << std::endl;
     }

This is better this way after all.

With this two modifications, you should be able to reproduc ESSL example. The output I got is:
Code: Select all
input A:
0       2
2       -1
2       -1
0       1.5
2       -1
2       -1

output A:
-4      2
0.5     2.5
0.5     0.285714
0       -0.428571
0.5     0.285714
0.5     0.285714


So that's the first point. Second point where are the factors Q and R.
First refer to the header of the LAPACK routine DGEQRF:
http://www.netlib.org/lapack/double/dgeqrf.f

R is the upper triangular part of A.
So for this example, we have
Code: Select all
R = [ 4 2
        0 2.5
        0 0
        0 0
        0 0
        0 0];


Q is more tricky to get. The lower part in A returns you the vector used in the Householder transformation, the vector tau returns you the coefficent for this transformation.

In our example:
Code: Select all
[ v_1  v_2 ]= [
                       1      0
                       0.5   1
                       0.5   0.285714
                       0    -0.428571
                       0.5   0.285714
                       0.5   0.285714 ];
[ tau_1 tau_2 ] = [ 1 1.4 ];

and
Code: Select all
H_1 = ( I - tau_1 *v_1 *v_1' )
H_2 = ( I - tau_2 *v_2 *v_2' )
Q = H_1 * H_2


You can check that with these formulas A= Q*R.

The routine DORMQR is quite handy routine if you want Q. You set the matrix C to identity in input and in output you will get Q if you set trans='N' and side='L'.

Julie.
Julie
 
Posts: 299
Joined: Wed Feb 23, 2005 12:32 am
Location: ICL, Denver. Colorado

lwork=0? dormqr

Postby wxiong » Mon Aug 28, 2006 6:15 pm

hi, julia:

Thanks for your help. It solves my puzzle.

I am using window’s version of Lapack. I thought it automatically changes to the row major from the column major during compile. Now I know it does not. I need to convert into column major before I call any Lapcak routine. But I only need to call twice before and after Lapack. For example, the following is right:

Input A in Row Major,
Conver A into Column Major,
Call DGEQRF
Call DORMQR
Call Ddgetrf
Call Dgetri
Conver A into Row Major.

But, the following is not correct:

Input A in Row Major,
Conver A into Column Major,
Call DGEQRF
Conver A into Column Major,
Call DORMQR
Conver A into Column Major,
Call Ddgetrf
Conver A into Column Major,
Call Dgetri
Conver A into Row Major.


Am I right?

Second question I have is about parameter lwork. In the link http://www.cepba.upc.edu/docs/ibm_doc/m ... sl368.html, it sayas “For best perfomance specify lwork = 0”. But, in my code, if I set lwork=0, I got compile error.

Although with your detailed explanations, I totally understood how to find or calculate Q and R, I like to know how to use routine DORMQR to find Q. I am not sure what is A and C. Is A the same input in DGEQRF’s A? is C the [V-1 V-2]?

Thanks for your generous help.

Wei.
wxiong
 
Posts: 3
Joined: Mon Aug 28, 2006 3:04 pm

Next

Return to User Discussion

Who is online

Users browsing this forum: No registered users and 3 guests

cron