CME for using both ATLAS & CLAPACK (zgesvd) in C++ ?

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

CME for using both ATLAS & CLAPACK (zgesvd) in C++ ?

Postby flexneck » Mon May 17, 2010 10:24 am

I am trying to integrate both ATLAS BLAS and CLAPACK routines together to implement an algorithm. I am having trouble getting the routines to link. My initial strategy was to create a routine using zgesvd, and a separate routine using ATLAS CBLAS functions. However, I'm having trouble getting the linker to bridge the difference between dcomp of type complex<double> and doublecomplex. I haven't found the correct recipe of where, or how to do the typecast. This is what I have so far...

/* svd.cpp */

#include <cstdlib>
#include <cmath>
#include <cstdio>
extern "C"{
#include "/home/bruce/Apps/CLAPACK-3.2.1/INCLUDE/f2c.h"
#include "/usr/local/atlas/include/cblas.h"
#include "/usr/local/atlas/include/clapack.h"
#include "/home/bruce/Apps/CLAPACK-3.2.1/INCLUDE/clapack.h"
#include "svd.h"

int mysvd( char jobu, char jobv, int m, int n, doublecomplex *A, int lda, double *S,
doublecomplex *U, int ldu, doublecomplex *VT, int ldvt, doublecomplex *wk, int lwork)
{
doublereal RWORK;
integer M = (integer)m;
integer N = (integer)n;
integer INFO;
integer LDA = (integer)lda;
integer LDU = (integer)ldu;
integer LDVT = (integer)ldvt;
integer LWORK = (integer)lwork;

zgesvd_( &jobu, &jobv, &M, &N, A, &LDA, S, U, &LDU, VT, &LDVT, wk, &LWORK, &RWORK, &INFO);
if ((long) INFO == 0 )
return 0;
else
return -1;
}
} /* end svd.cpp */

/* Function Prototype for mysvd, svd.h */
#ifdef __cplusplus
extern "C"{
#endif
//#include "/home/bruce/Apps/CLAPACK-3.2.1/INCLUDE/f2c.h"

#ifndef doublecomplex
typedef struct{ double r; double i;} doublecomplex;
#endif

int mysvd( char jobu, char jobv, int m, int n, doublecomplex *A, int lda, double *S,
doublecomplex *U, int ldu, doublecomplex *VT, int ldvt,
doublecomplex *wk, int lwork);

#ifdef __cplusplus
}
#endif

/*
** test.cpp -- an initial stab at integrating CBLAS & CLAPACK functions
This routine uses both BLAS and CLAPACK libraries from netlib.org
*/
/* includes go here */
#include <iostream>
#include <algorithm>
#include <stdio.h>
#include <sys/time.h>
#include <time.h>
#include <complex>

/* put in my includes here */

using namespace std;
using std::complex;
typedef complex<double> dcomp; /* Define complex double data type */

extern "C"{ /* this avoids name mangling */
#include "/usr/local/atlas/include/cblas.h"
#include "/usr/local/atlas/include/clapack.h"
#include "svd.h"
}

void MAIN_(){}
void MAIN__(){}
void _MAIN_(){}

int main(void)
{
long m,n,ii,jj;
double diff, NN;
double realsig, imagsig;

m = 2; n = 2;
//int mm,kk,nn; // matrix dimensions, A = [mm x kk], B = [kk x nn]
//mm = m; kk = n; nn = m;
int ldA, ldB, ldC;
ldA = m; ldB = m; ldC = m;
dcomp A[m][m];
dcomp C[m][m];
dcomp I[m][m];
/* make identity matrix */
for ( ii=0; ii<m; ii++){
for (jj=0; jj<m; jj++){
if (ii == jj)
I[ii][jj] = dcomp(1.0, 0.0);
else
I[ii][jj] = dcomp(0.0, 0.0);
}
}

/* A is a matrix which is the result of a bunch of CBLAS routines. To
simplify, let us start with a matrix initialize to the following
*/

/* below is Fortran Order. Previous routines used C order, so we will need
to convert for use in svd
A[0].r = 0.69096145; A[0].i = 0.54390208; // A[0][0]
A[1].r = 0.83805603; A[1].i = 0.27186403; // A[1][0]
A[2].r = 0.35559057; A[2].i = 0.12026953; // A[0][1]
A[3].r = 0.89833702; A[3].i = 0.06791132; // A[1][1]
*/
/* Assume routine starts with CBLASRowMajor */
A[0][0] = dcomp( 0.69096145, 0.54390208 ); // A[0][0] A[0]
A[0][1] = dcomp( 0.35559057, 0.12026953 ); // A[0][1] A[1]
A[1][0] = dcomp( 0.83805603, 0.27186403 ); // A[1][0] A[2]
A[1][1] = dcomp( 0.89833702, 0.06791132 ); // A[1][1] A[3]

/* Now form the transpose using zgemm */
dcomp alpha = dcomp(1.0, 0.0);
dcomp beta = dcomp(0.0, 0.0);

cblas_zgemm( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, m,
&alpha, A, ldA, I, ldB, &beta, C, ldC);

/* If all goes well, C will now contain A transpose */

/* print it out to check... */
for (ii=0; ii<m; ii++){
for (jj=0; jj<m; jj++){
printf("C[%li][%li]= %12.8lg, %12.8lgj\n", ii, jj,
real(C[ii][jj]), imag(C[ii][jj]));
}
}

/* define variables for mysvd */

char JOBU, JOBVT;
JOBU = 'A'; JOBVT = 'A';
int lda = m; int ldu = m;
int ldvt = n;

double S[m];
dcomp UU[m][m];
dcomp VV[m][m];

int mn = min(m,n);
int MN = max(m,n);
int lwork = 2*max(3*mn+MN,5*mn);
dcomp wk[lwork];
//typedef struct{ double r; double i;} doublecomplex;

//if ( mysvd( JOBU, JOBVT, m, n, (doublecomplex*)&C[0][0], lda, &S[0], (doublecomplex*)&UU[0][0], ldu, (doublecomplex*)&VV[0][0], ldvt, (doublecomplex*)&wk[0], lwork) == 0)
if ( mysvd( JOBU, JOBVT, m, n, C, lda, S, UU, ldu, VV, ldvt, wk, lwork) == 0)
//if ( mysvd( JOBU, JOBVT, m, n, *C, lda, *S, *UU, ldu, *VV, ldvt, *wk, lwork) == 0)
{
cout << "Hey that seemed to work! Let's continue development..." << endl;
}
else{
cout << "bummer, this attempt did not work" << endl;
}

return 0;
}

Makefile for above

# Makefile for test

CC = g++

CFLAGS = -O3 -m64 -Wall -Wcast-qual

ROOTPATH = /home/bruce/Apps/CLAPACK-3.2.1

BLASPATH = /usr/local/atlas

INCDIRS = -I$(BLASPATH)/include/ -I$(ROOTPATH)/INCLUDE/ -I$(ROOTPATH)

INCSVDDIRS = $(INCDIRS)

F2CDIR = $(ROOTPATH)/F2CLIBS

LDLIBS = $(ROOTPATH)/lapack_LINUX.a $(ROOTPATH)/libcblaswr.a $(BLASPATH)/lib/libcblas.a $(BLASPATH)/lib/libatlas.a $(F2CDIR)/libf2c.a -lm

test: test.o svd.o
$(CC) test.o svd.o $(LDLIBS) -o esprit

svd.o: svd.cpp
$(CC) svd.cpp $(INCSVDDIRS) -c

test.o: test.cpp
$(CC) test.cpp $(INCDIRS) -c

clean:
rm -f test test.o svd.o

$ make test
g++ test.cpp -I/usr/local/atlas/include/ -I/home/bruce/Apps/CLAPACK-3.2.1/INCLUDE/ -I/home/bruce/Apps/CLAPACK-3.2.1 -c
test.cpp: In function ‘int main()’:
test.cpp:105: error: cannot convert ‘dcomp (*)[(((long unsigned int)(((long int)m) + -0x00000000000000001)) + 1)]’ to ‘doublecomplex*’ for argument ‘5’ to ‘int mysvd(char, char, int, int, doublecomplex*, int, double*, doublecomplex*, int, doublecomplex*, int, doublecomplex*, int)’
make: *** [test.o] Error 1

Various permutations of modifying the files have not resulted in anything much better. Can someone offer some insight on how to bridge the C++ world over to the f2c world? I had a very experienced C programmer sit with me last week and try literally dozens of ideas, none of which panned out.

I thought I'd post this to the list to get some tips, does, & don'ts on this sort of thing. Thanks for any insights!
flexneck
 
Posts: 7
Joined: Thu May 06, 2010 10:25 am

Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot] and 2 guests