Backing up a little - Trying to get LAPACK to work...

bruce.labitt at autoliv.com bruce.labitt at autoliv.com
Wed May 26 14:19:33 EDT 2010


gnhlug-discuss-bounces at mail.gnhlug.org wrote on 05/26/2010 09:10:29 AM:

> 
> How are you invoking Valgrind?  Where is the Valgrind output for the
> 9x9 run?
> 
> --kevin
> -- 
> alumni.unh.edu!kdc / http://kdc-blog.blogspot.com/
> GnuPG: D87F DAD6 0291 289C EB1E 781C 9BF8 A7D8 B280 F24E
> 
>  Wipe him down with gasoline 'til his arms are hard and mean
>  From now on boys this iron boat's your home
>  So heave away, boys.
>    -- Tom Waits
> 
> _______________________________________________
>

here you go...
last post to the list got bounced - too big, now with some edits to keep 
char count down
======= main program ========

/* test1.c, a simple shell program to call svd.c */
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <time.h>
#include <ctype.h>
#include <stdarg.h>
#include <string.h>
#include <stdint.h>
#include <inttypes.h>
#include <math.h>

/* CBLAS was created by ATLAS using gcc and gfortran */
#include "/usr/local/atlas/include/cblas.h"
#include "/usr/local/atlas/include/clapack.h"
#include "svd.h"

//#define littletry
//#define onstack
//#define mylwork
#define moveme
#define setrlim
/* ifdef littletry, use a 2 x 2 Complex Matrix */

/* some FORTRAN stuff */
void MAIN_(){};
void MAIN__(){};
void _MAIN_(){};


int main(void)
{ 
    long m,n,ii,jj;
    //double diff, NN;
    //double realsig, imagsig;
    #ifdef littletry
        m = 2; n = 2;
        #else
        m = 9; n = 9; 
        #endif 
 
    #ifdef setrlim
    const rlim_t kStackSize = 2 * 1024 * 1024;   // min stack size = 16 MB
    struct rlimit rl;
    int result;

    result = getrlimit(RLIMIT_STACK, &rl);
    if (result == 0)
    {
        if (rl.rlim_cur < kStackSize)
        {
            rl.rlim_cur = kStackSize;
            result = setrlimit(RLIMIT_STACK, &rl);
            if (result != 0)
            {
                fprintf(stderr, "setrlimit returned result = %d\n", 
result);
            }
        }
    }
    #endif
 
    #ifdef moveme
    double S[m];                                        /* in MATLAB, this 
is D */
    for (ii=0; ii<m; ii++){
        S[ii] = 0.0;
    }
    #ifdef onstack
    doublecomplex UU[m*m];
    doublecomplex VV[m*m];
    #else
    doublecomplex *UU;  doublecomplex *VV;  doublecomplex *wk;
    UU = (doublecomplex*)malloc(m*m*sizeof(doublecomplex));
    VV = (doublecomplex*)malloc(m*m*sizeof(doublecomplex));
    for (ii=0; ii<m*m; ii++){
        UU[ii].r = 0.0; UU[ii].i = 0.0;
    }
    for (ii=0; ii<n*n; ii++){
        VV[ii].r = 0.0; VV[ii].i = 0.0;
    }
    #endif
    #endif
 
    int ldA, ldB, ldC;
    ldA = m; ldB = m; ldC = m;
 
    #ifdef onstack
    doublecomplex A[m*n];  doublecomplex C[m*n];  doublecomplex I[m*n];
    #else
    doublecomplex *A;  doublecomplex *C;  doublecomplex *I;
    A = (doublecomplex*) malloc( m * n * sizeof(doublecomplex));
    C = (doublecomplex*) malloc( m * n * sizeof(doublecomplex));
    I = (doublecomplex*) malloc( m * n * sizeof(doublecomplex));
    #endif
 
    /* make identity matrix */
    for ( ii=0; ii<m; ii++){
        for (jj=0; jj<m; jj++){
            if (ii == jj){
                I[ii*n+jj].r = 1.0; I[ii*n+jj].i = 0.0;}
            else{
                I[ii*n+jj].r = 0.0; I[ii*n+jj].i = 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 */
 
    /* original data, which works */
    #ifdef littletry
    /* a 2 x 2 */
    A[0*n + 0].r = 0.69096145; A[0*n + 0].i = 0.54390208; // A[0][0]  A[0]
    A[0*n + 1].r = 0.35559057; A[0*n + 1].i = 0.12026953; // A[0][1]  A[1]
    A[1*n + 0].r = 0.83805603; A[1*n + 0].i = 0.27186403; // A[1][0]  A[2]
    A[1*n + 1].r = 0.89833702; A[1*n + 1].i = 0.06791132; // A[1][1]  A[3]
    #else
    /* a 9 x 9 */
    A[0*n + 0].r =    0.1056245;  A[0*n + 0].i = 3.008661e-18;
        A[0*n + 1].r =   0.12341259;  A[0*n + 1].i = -0.048803565;
        A[0*n + 2].r =   0.12811331;  A[0*n + 2].i =  -0.10658785;
        A[0*n + 3].r =   0.10972881;  A[0*n + 3].i =  -0.16601482;
        A[0*n + 4].r =  0.066265673;  A[0*n + 4].i =  -0.21196376;
        A[0*n + 5].r = 0.0050456614;  A[0*n + 5].i =  -0.23662813;
        A[0*n + 6].r = -0.060823528;  A[0*n + 6].i =  -0.23785819;
        A[0*n + 7].r =  -0.12032375;  A[0*n + 7].i =  -0.22090523;
        A[0*n + 8].r =  -0.17269674;  A[0*n + 8].i =  -0.18375118;
        A[1*n + 0].r =   0.12341259;  A[1*n + 0].i =  0.048803565;
        A[1*n + 1].r =   0.23163068;  A[1*n + 1].i = 1.6263033e-17;
        A[1*n + 2].r =   0.33062898;  A[1*n + 2].i = -0.091557204;
        A[1*n + 3].r =   0.39583681;  A[1*n + 3].i =  -0.22168363;
        A[1*n + 4].r =   0.40637821;  A[1*n + 4].i =  -0.36824733;
        A[1*n + 5].r =   0.36179959;  A[1*n + 5].i =  -0.51181571;
        A[1*n + 6].r =   0.26657764;  A[1*n + 6].i =  -0.63062745;
        A[1*n + 7].r =   0.13972331;  A[1*n + 7].i =  -0.71466428;
        A[1*n + 8].r = -0.011894371;  A[1*n + 8].i =  -0.75281686;
        A[2*n + 0].r =   0.12811331;  A[2*n + 0].i =   0.10658785;
        A[2*n + 1].r =   0.33062898;  A[2*n + 1].i =  0.091557204;
        A[2*n + 2].r =   0.54440918;  A[2*n + 2].i = 6.6542908e-18;
        A[2*n + 3].r =    0.7242013;  A[2*n + 3].i =  -0.17340506;
        A[2*n + 4].r =    0.8296598;  A[2*n + 4].i =  -0.40593058;
        A[2*n + 5].r =   0.84813909;  A[2*n + 5].i =  -0.66878153;
        A[2*n + 6].r =   0.76642337;  A[2*n + 6].i =  -0.92313234;
        A[2*n + 7].r =   0.60716441;  A[2*n + 7].i =   -1.1398751;
        A[2*n + 8].r =   0.37907731;  A[2*n + 8].i =   -1.2956598;
        A[3*n + 0].r =   0.10972881;  A[3*n + 0].i =   0.16601482;
        A[3*n + 1].r =   0.39583681;  A[3*n + 1].i =   0.22168363;
        A[3*n + 2].r =    0.7242013;  A[3*n + 2].i =   0.17340506;
        A[3*n + 3].r =    1.0389688;  A[3*n + 3].i = 1.8458542e-17;
        A[3*n + 4].r =    1.2733233;  A[3*n + 4].i =  -0.28371476;
        A[3*n + 5].r =    1.4011977;  A[3*n + 5].i =  -0.64373949;
        A[3*n + 6].r =    1.3852911;  A[3*n + 6].i =   -1.0307846;
        A[3*n + 7].r =     1.243208;  A[3*n + 7].i =   -1.3947563;
        A[3*n + 8].r =    0.9826123;  A[3*n + 8].i =   -1.6993985;
        A[4*n + 0].r =  0.066265673;  A[4*n + 0].i =   0.21196376;
        A[4*n + 1].r =   0.40637821;  A[4*n + 1].i =   0.36824733;
        A[4*n + 2].r =    0.8296598;  A[4*n + 2].i =   0.40593058;
        A[4*n + 3].r =    1.2733233;  A[4*n + 3].i =   0.28371476;
        A[4*n + 4].r =    1.6548532;  A[4*n + 4].i = 1.9515639e-18;
        A[4*n + 5].r =     1.925963;  A[4*n + 5].i =  -0.41275362;
        A[4*n + 6].r =    2.0247087;  A[4*n + 6].i =  -0.90354436;
        A[4*n + 7].r =    1.9556722;  A[4*n + 7].i =   -1.4043293;
        A[4*n + 8].r =    1.7198514;  A[4*n + 8].i =   -1.8659826;
        A[5*n + 0].r = 0.0050456614;  A[5*n + 0].i =   0.23662813;
        A[5*n + 1].r =   0.36179959;  A[5*n + 1].i =   0.51181571;
        A[5*n + 2].r =   0.84813909;  A[5*n + 2].i =   0.66878153;
        A[5*n + 3].r =    1.4011977;  A[5*n + 3].i =   0.64373949;
        A[5*n + 4].r =     1.925963;  A[5*n + 4].i =   0.41275362;
        A[5*n + 5].r =    2.3602752;  A[5*n + 5].i = 1.5612511e-17;
        A[5*n + 6].r =    2.6099806;  A[5*n + 6].i =  -0.55240669;
        A[5*n + 7].r =    2.6629588;  A[5*n + 7].i =   -1.1625421;
        A[5*n + 8].r =    2.5083167;  A[5*n + 8].i =   -1.7713606;
        A[6*n + 0].r = -0.060823528;  A[6*n + 0].i =   0.23785819;
        A[6*n + 1].r =   0.26657764;  A[6*n + 1].i =   0.63062745;
        A[6*n + 2].r =   0.76642337;  A[6*n + 2].i =   0.92313234;
        A[6*n + 3].r =    1.3852911;  A[6*n + 3].i =    1.0307846;
        A[6*n + 4].r =    2.0247087;  A[6*n + 4].i =   0.90354436;
        A[6*n + 5].r =    2.6099806;  A[6*n + 5].i =   0.55240669;
        A[6*n + 6].r =    3.0257882;  A[6*n + 6].i = -1.9125326e-16;
        A[6*n + 7].r =    3.2336681;  A[6*n + 7].i =  -0.66610971;
        A[6*n + 8].r =    3.2113302;  A[6*n + 8].i =   -1.3815207;
        A[7*n + 0].r =  -0.12032375;  A[7*n + 0].i =   0.22090523;
        A[7*n + 1].r =   0.13972331;  A[7*n + 1].i =   0.71466428;
        A[7*n + 2].r =   0.60716441;  A[7*n + 2].i =    1.1398751;
        A[7*n + 3].r =     1.243208;  A[7*n + 3].i =    1.3947563;
        A[7*n + 4].r =    1.9556722;  A[7*n + 4].i =    1.4043293;
        A[7*n + 5].r =    2.6629588;  A[7*n + 5].i =    1.1625421;
        A[7*n + 6].r =    3.2336681;  A[7*n + 6].i =   0.66610971;
        A[7*n + 7].r =    3.6083259;  A[7*n + 7].i = 1.8648277e-16;
        A[7*n + 8].r =    3.7462007;  A[7*n + 8].i =  -0.77126832;
        A[8*n + 0].r =  -0.17269674;  A[8*n + 0].i =   0.18375118;
        A[8*n + 1].r = -0.011894371;  A[8*n + 1].i =   0.75281686;
        A[8*n + 2].r =   0.37907731;  A[8*n + 2].i =    1.2956598;
        A[8*n + 3].r =    0.9826123;  A[8*n + 3].i =    1.6993985;
        A[8*n + 4].r =    1.7198514;  A[8*n + 4].i =    1.8659826;
        A[8*n + 5].r =    2.5083167;  A[8*n + 5].i =    1.7713606;
        A[8*n + 6].r =    3.2113302;  A[8*n + 6].i =    1.3815207;
        A[8*n + 7].r =    3.7462007;  A[8*n + 7].i =   0.77126832;
        A[8*n + 8].r =    4.0597545;  A[8*n + 8].i = -1.1709383e-17;
        #endif
 
        printf("\nOriginal matrix in CblasRowMajor form\n\n");
        for (ii=0; ii<m; ii++){
        for (jj=0; jj<n; jj++){
            printf("A[%li][%li]==> A[%li] = %12.8lg, %12.8lgj\n", ii, jj, 
                     ii*n + jj, A[ii*n + jj].r, A[ii*n + jj].i);
        }
    }
        /* Now form the transpose using zgemm */
        doublecomplex alpha; 
        alpha.r = 1.000000; alpha.i = 0.000000;
        doublecomplex beta; 
        beta.r  = 0.000000; beta.i  = 0.000000;
 
        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... */
    printf("\nNow matrix A is transposed by zgemm, but the indices are 
strange\n\n");
    for (ii=0; ii<m; ii++){
        for (jj=0; jj<n; jj++){
            printf("C[%li][%li]==> C[%li] = %12.8lg, %12.8lgj\n", ii, jj, 
                 ii*n+jj, C[ii*n + jj].r, C[ii*n + jj].i);
        }
    }
 
    printf("\nHowever, to LAPACK it looks like:\n");
    for (jj=0; jj<n; jj++){
        for (ii=0; ii<m; ii++){
            printf("C[%li][%li]==> C[%li] = %12.8lg, %12.8lgj\n", ii, jj, 
                       ii + m * jj, C[ii + m * jj].r, C[ii+ m*jj].i);
        }
    }
    printf("\n");
 
    #ifndef onstack
    free(A);
    free(I);
    #endif
 
        /* define variables for mysvd */
 
        char JOBU, JOBVT;
        JOBU = 'A';     JOBVT = 'A';
        int lda = m;    int ldu = m;
        int ldvt = n;
 
        #ifndef moveme
        double S[m];                                    /* in MATLAB, this 
is D */
        for (ii=0; ii<m; ii++){
        S[ii] = 0.0;
    }
        #ifdef onstack
        doublecomplex UU[m*m];
        doublecomplex VV[m*m];
        #else
        doublecomplex *UU;  doublecomplex *VV;  doublecomplex *wk;
        UU = (doublecomplex*)malloc(m*m*sizeof(doublecomplex));
        VV = (doublecomplex*)malloc(m*m*sizeof(doublecomplex));
    for (ii=0; ii<m*m; ii++){
        UU[ii].r = 0.0; UU[ii].i = 0.0;
    }
    for (ii=0; ii<n*n; ii++){
        VV[ii].r = 0.0; VV[ii].i = 0.0;
    } 
        #endif
        #endif
 
        int mn = min(m,n); 
        int MN = max(m,n);

        int lwork = 2*max(1, 2*mn+MN );     // this is the wrong size...
        #ifdef mylwork
        lwork = 891;
        #endif
        #ifdef onstack
        doublecomplex wk[lwork];
        #else
        wk = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
    #endif
 
        int good = mysvd( JOBU, JOBVT, m, n, &C[0], lda, &S[0], &UU[0], 
ldu, &VV[0], 
                                                        ldvt, &wk[0], 
lwork);
        if (good == 0) { 
                printf("\nHey that seemed to work!  Let's continue 
development...\n");
                printf("\nThe matrix UU is:\n");
                for (jj=0;jj<n; jj++){
                        for (ii=0; ii<m; ii++){
                                printf("UU[%li][%li]==> UU[%li] = %12.8lg, 
%12.8lgj\n", ii, jj, 
                       ii+m*jj, UU[ii*n + jj].r, UU[ii*n + jj].i);
                }
        }
                printf("\nThe matrix V is\n");
                for (jj=0;jj<n; jj++){
                        for (ii=0; ii<m; ii++){
                                printf("V[%li][%li]==> VV[%li] = %12.8lg, 
%12.8lgj\n", ii, jj, 
                       ii+m*jj, VV[ii*n + jj].r, VV[ii*n + jj].i);
                }
        }
        printf("\nThe array S is\n");
        for (ii=0; ii<m; ii++){
                        printf("S[%li]= %12.8lg\n", ii, S[ii]);
        }
        }
        else { printf("bummer, this attempt did not work\n"); }
        #ifndef onstack
        free(UU); free(VV); free(C); free(wk);
        #endif
        return 0;
}

====== svd.c  actually calls zgesvd_ 
======================================

/* svd.c */

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

#include "f2c.h"
#include "/usr/local/atlas/include/cblas.h"
#include "/usr/local/atlas/include/clapack.h"

//#define mydebug
#define real_lwork  // means we we sent a legitimate value
#define LAPACK      // linked with LAPACK, not CLAPACK

/* some FORTRAN flooby dust needed for CLAPACK */
#ifndef LAPACK
void MAIN_(){}
void MAIN__(){}
void _MAIN_(){}
#endif

/* the following is the convention to call an external FORTRAN program */

extern void (zgesvd_)( char *jobu, char *jobvt, long *m, long *n, 
doublecomplex *A, 
    long *lda, double *s, doublecomplex *u, long *ldu, doublecomplex *vt, 
long *ldvt, 
    doublecomplex *work, long *lwork, double *rwork, long *info);

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;
        long myinfo = -32;
    int retval;
    //long ii, jj;
 
    printf("Got inside svd.c, just before zgesvd_ call.\n");
    printf("jobu = %c, jobv = %c, m = %li, n = %li, LDA = %li, LDU = %li, 
\
    LDVT = %li, LWORK = %li\n", jobu, jobv, (long)m, (long)n, (long)LDA, \
    (long)LDU, (long)LDVT, (long)LWORK);
    /* 
    zgesvd_( &jobu, &jobv, &M, &N, A, &LDA, S, U, &LDU, VT, &LDVT, wk, 
                 &LWORK, &RWORK, &INFO);
    */
 
    LWORK = -1;
    zgesvd_( &jobu, &jobv, &M, &N, A, &LDA, S, U, &LDU, VT, &LDVT, wk, 
                 &LWORK, &RWORK, &INFO);
 
    /* call svd to figure out how big LWORK really should be */
    myinfo = (long)INFO & 0x00000000FFFFFFFF;
    printf("myinfo = %li\n", myinfo);
 
    #ifdef mydebug
    /* write to U, V, A */
    for (ii=0; ii< m; ii++){
        for (jj=0; jj<n; jj++){
            U[ii*n + jj].r  =  1.0000002; U[ii*n + jj].i = -1.0001;
            VT[ii*n + jj].r = -1.01; VT[ii*n + jj].i  =  1.0000000003;
        }
    }
    for (ii=0; ii< m; ii++){
        for (jj=0; jj<n; jj++){
            A[ii*n + jj].r = 0.1;  A[ii*n + jj].i = -0.6;
        }
    }
    for (ii=0; ii<m; ii++){
        S[ii] = m-ii;
    }
    printf("optimal LWORK = wk[0].r = %lg\n", wk[0].r);
    printf("now checking wk[0].i too = %lg\n", wk[0].i);
    LWORK = (integer)(wk[0].r*1.1);
    printf("New LWORK = wk[0].r * 1.1 = %li\n", (long)(1.1*wk[0].r));
    doublecomplex *WK;
    WK = (doublecomplex*) malloc ( sizeof(doublecomplex) * LWORK );
    printf("Allocated to WK, LWORK = %li values\n", (long)(1.1*wk[0].r));
    for (ii=0; ii< (long)(1.1*wk[0].r); ii++){
        WK[ii].r = ii; WK[ii].i = ii-3.14159;
        //printf("WK[%li] = %lg + %lgj\n", ii, WK[ii].r, WK[ii].i);
    }
    free(WK);
    printf("WK freed\n");
    retval = 0;
    return retval;
 
    #else                       // we are NOT in mydebug.  We are LIVE!
    if (myinfo == 0){
        printf("optimal LWORK = wk[0].r = %lg\n", wk[0].r);
        printf("now checking wk[0].i too = %lg\n", wk[0].i);
        LWORK = (integer)(wk[0].r*1.2);
        printf("New LWORK = wk[0].r * 1.2 = %li\n", (long)LWORK);
        doublecomplex *WK, *junk;
        WK   = (doublecomplex*) malloc ( sizeof(doublecomplex) * LWORK );
        junk = (doublecomplex*) malloc ( sizeof(doublecomplex) * LWORK * 
10000 );
        printf("Just before the second call to zgesvd\n");
        zgesvd_( &jobu, &jobv, &M, &N, A, &LDA, S, U, &LDU, VT, &LDVT, WK, 

                 &LWORK, &RWORK, &INFO);
        printf("Just after the second call to zgesvd\n");
        myinfo = (long)INFO & 0x00000000FFFFFFFF;
        printf("INFO = %li\n", myinfo);
        free(junk);
        printf("junk was freed\n");
        /*
        for (ii=0; ii<m; ii++){
            printf("S[%li] = %lg\n", ii, S[ii]);
            }
            for (jj=0; jj<n; jj++){
                for (ii=0; ii<m; ii++){
                    printf("U[%li][%li] ==> U[%li] = %lg + %lgj\n", ii, 
jj, ii+m*jj, U[ii+m*jj].r, U[ii+m*jj].i);
            }
        }
        */
        free(WK);
        retval = 0;
        return retval;
        }
    else
    {
        printf("First call with LWORK=-1 returned with an error\n");
        retval = -1;
        return retval;
        }
        /*
    if (myinfo == 0 )
    {
        printf("The SVD was successful, now exiting svd.c\n");
        retval = 0;
            return retval;
    }
    else
    { 
        printf("Result of second call has some kind of error\n");
        retval = -1;
        return retval; 
    }
    */
    #endif
}

========== svd.h ============================================

/* Function Prototype for mysvd, svd.h */

//#ifdef __cplusplus
//extern "C"{
//#endif
#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

========= Makefile ===========================================

# Makefile for test1 with LAPACK

CC = gcc

CFLAGST = -O3 -m64 -mtune=native -march=native -g -v -fstack-protector-all 
-Wall -Wcast-qual
CFLAGSS = -O3 -m64 -mtune=native -march=native -g -v -fstack-protector-all 
-Wall -Wcast-qual

LINKERFLAGS = -Wl,-Map=test1.map,--cref

ROOTPATH = /home/bruce/APPS/lapack-3.2.1

BLASPATH = /usr/local/atlas

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

INCSVDDIRS = $(INCDIRS) -I.

#F2CDIR  = $(ROOTPATH)/F2CLIBS

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

LDLIBS = $(ROOTPATH)/lapack_LINUX.a $(BLASPATH)/lib/libcblas.a 
$(BLASPATH)/lib/libf77blas.a $(BLASPATH)/lib/liblapack.a 
$(BLASPATH)/lib/libatlas.a -lgfortran -lm

test1: test1.o svd.o
        $(CC) $(LINKERFLAGS) test1.o svd.o $(LDLIBS) -o test1

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

test1.o:        test1.c
        $(CC) $(CFLAGSS) test1.c $(INCDIRS) -c
 
clean: 
        rm -f test1 test1.o svd.o

========== And finally Valgrind Output 
========================================================

$ valgrind --leak-check=full --show-reachable=yes ./test1
==2524== Memcheck, a memory error detector
==2524== Copyright (C) 2002-2009, and GNU GPL'd, by Julian Seward et al.
==2524== Using Valgrind-3.6.0.SVN-Debian and LibVEX; rerun with -h for 
copyright info
==2524== Command: ./test1
==2524== 

Original matrix in CblasRowMajor form

A[0][0]==> A[0] =    0.1056245, 3.008661e-18j
A[0][1]==> A[1] =   0.12341259, -0.048803565j
A[0][2]==> A[2] =   0.12811331,  -0.10658785j
A[0][3]==> A[3] =   0.10972881,  -0.16601482j
A[0][4]==> A[4] =  0.066265673,  -0.21196376j
A[0][5]==> A[5] = 0.0050456614,  -0.23662813j
A[0][6]==> A[6] = -0.060823528,  -0.23785819j
A[0][7]==> A[7] =  -0.12032375,  -0.22090523j
A[0][8]==> A[8] =  -0.17269674,  -0.18375118j
A[1][0]==> A[9] =   0.12341259,  0.048803565j
A[1][1]==> A[10] =   0.23163068, 1.6263033e-17j
A[1][2]==> A[11] =   0.33062898, -0.091557204j
A[1][3]==> A[12] =   0.39583681,  -0.22168363j
A[1][4]==> A[13] =   0.40637821,  -0.36824733j
A[1][5]==> A[14] =   0.36179959,  -0.51181571j
A[1][6]==> A[15] =   0.26657764,  -0.63062745j
A[1][7]==> A[16] =   0.13972331,  -0.71466428j
A[1][8]==> A[17] = -0.011894371,  -0.75281686j
A[2][0]==> A[18] =   0.12811331,   0.10658785j
A[2][1]==> A[19] =   0.33062898,  0.091557204j
A[2][2]==> A[20] =   0.54440918, 6.6542908e-18j
A[2][3]==> A[21] =    0.7242013,  -0.17340506j
A[2][4]==> A[22] =    0.8296598,  -0.40593058j
A[2][5]==> A[23] =   0.84813909,  -0.66878153j
A[2][6]==> A[24] =   0.76642337,  -0.92313234j
A[2][7]==> A[25] =   0.60716441,   -1.1398751j
A[2][8]==> A[26] =   0.37907731,   -1.2956598j
A[3][0]==> A[27] =   0.10972881,   0.16601482j
A[3][1]==> A[28] =   0.39583681,   0.22168363j
A[3][2]==> A[29] =    0.7242013,   0.17340506j
A[3][3]==> A[30] =    1.0389688, 1.8458542e-17j
A[3][4]==> A[31] =    1.2733233,  -0.28371476j
A[3][5]==> A[32] =    1.4011977,  -0.64373949j
A[3][6]==> A[33] =    1.3852911,   -1.0307846j
A[3][7]==> A[34] =     1.243208,   -1.3947563j
A[3][8]==> A[35] =    0.9826123,   -1.6993985j
<snip>


However, to LAPACK it looks like:
C[0][0]==> C[0] =    0.1056245, 3.008661e-18j
C[1][0]==> C[1] =   0.12341259,  0.048803565j
C[2][0]==> C[2] =   0.12811331,   0.10658785j
C[3][0]==> C[3] =   0.10972881,   0.16601482j
C[4][0]==> C[4] =  0.066265673,   0.21196376j
C[5][0]==> C[5] = 0.0050456614,   0.23662813j
C[6][0]==> C[6] = -0.060823528,   0.23785819j
C[7][0]==> C[7] =  -0.12032375,   0.22090523j
C[8][0]==> C[8] =  -0.17269674,   0.18375118j
C[0][1]==> C[9] =   0.12341259, -0.048803565j
C[1][1]==> C[10] =   0.23163068, 1.6263033e-17j
C[2][1]==> C[11] =   0.33062898,  0.091557204j
C[3][1]==> C[12] =   0.39583681,   0.22168363j
C[4][1]==> C[13] =   0.40637821,   0.36824733j
C[5][1]==> C[14] =   0.36179959,   0.51181571j
C[6][1]==> C[15] =   0.26657764,   0.63062745j
C[7][1]==> C[16] =   0.13972331,   0.71466428j
C[8][1]==> C[17] = -0.011894371,   0.75281686j
C[0][2]==> C[18] =   0.12811331,  -0.10658785j
C[1][2]==> C[19] =   0.33062898, -0.091557204j
C[2][2]==> C[20] =   0.54440918, 6.6542908e-18j
C[3][2]==> C[21] =    0.7242013,   0.17340506j
C[4][2]==> C[22] =    0.8296598,   0.40593058j
C[5][2]==> C[23] =   0.84813909,   0.66878153j
C[6][2]==> C[24] =   0.76642337,   0.92313234j
C[7][2]==> C[25] =   0.60716441,    1.1398751j
C[8][2]==> C[26] =   0.37907731,    1.2956598j
C[0][3]==> C[27] =   0.10972881,  -0.16601482j
C[1][3]==> C[28] =   0.39583681,  -0.22168363j
C[2][3]==> C[29] =    0.7242013,  -0.17340506j
C[3][3]==> C[30] =    1.0389688, 1.8458542e-17j
C[4][3]==> C[31] =    1.2733233,   0.28371476j
C[5][3]==> C[32] =    1.4011977,   0.64373949j
C[6][3]==> C[33] =    1.3852911,    1.0307846j
C[7][3]==> C[34] =     1.243208,    1.3947563j
C[8][3]==> C[35] =    0.9826123,    1.6993985j
C[0][4]==> C[36] =  0.066265673,  -0.21196376j
C[1][4]==> C[37] =   0.40637821,  -0.36824733j
C[2][4]==> C[38] =    0.8296598,  -0.40593058j
C[3][4]==> C[39] =    1.2733233,  -0.28371476j
C[4][4]==> C[40] =    1.6548532, 1.9515639e-18j
C[5][4]==> C[41] =     1.925963,   0.41275362j
C[6][4]==> C[42] =    2.0247087,   0.90354436j
C[7][4]==> C[43] =    1.9556722,    1.4043293j
C[8][4]==> C[44] =    1.7198514,    1.8659826j
<snip>

Got inside svd.c, just before zgesvd_ call.
jobu = A, jobv = A, m = 9, n = 9, LDA = 9, LDU = 9,     LDVT = 9, LWORK = 
54
myinfo = 0
optimal LWORK = wk[0].r = 594
now checking wk[0].i too = 0
New LWORK = wk[0].r * 1.2 = 712
Just before the second call to zgesvd
Just after the second call to zgesvd
INFO = 0
junk was freed

Hey that seemed to work!  Let's continue development...

The matrix UU is:
==2524== Invalid read of size 8
==2524==    at 0x408C7F: main (in 
/home/bruce/sb/Private/Bruce_Labitt/ESPRIT_prod/C/esprit_proj/test_integration/test1)
==2524==  Address 0x8 is not stack'd, malloc'd or (recently) free'd
==2524== 
==2524== 
==2524== Process terminating with default action of signal 11 (SIGSEGV)
==2524==  Access not within mapped region at address 0x8
==2524==    at 0x408C7F: main (in 
/home/bruce/sb/Private/Bruce_Labitt/ESPRIT_prod/C/esprit_proj/test_integration/test1)
==2524==  If you believe this happened as a result of a stack
==2524==  overflow in your program's main thread (unlikely but
==2524==  possible), you can try to increase the size of the
==2524==  main thread stack using the --main-stacksize= flag.
==2524==  The main thread stack size used in this run was 8388608.
==2524== 
==2524== HEAP SUMMARY:
==2524==     in use at exit: 7,321 bytes in 19 blocks
==2524==   total heap usage: 69 allocs, 50 frees, 113,946,345 bytes 
allocated
==2524== 
==2524== 5 bytes in 1 blocks are still reachable in loss record 1 of 19
==2524==    at 0x4C274A8: malloc (vg_replace_malloc.c:236)
==2524==    by 0x4E453D8: ??? (in /usr/lib/libgfortran.so.3.0.0)
==2524==    by 0x4EED9EB: ??? (in /usr/lib/libgfortran.so.3.0.0)
==2524==    by 0x4E45287: ??? (in /usr/lib/libgfortran.so.3.0.0)
==2524==    by 0x4F01725: ??? (in /usr/lib/libgfortran.so.3.0.0)
==2524==    by 0x4E41D5A: ??? (in /usr/lib/libgfortran.so.3.0.0)
==2524== 

<snip>

==2524== 
==2524== 864 bytes in 1 blocks are definitely lost in loss record 16 of 19
==2524==    at 0x4C274A8: malloc (vg_replace_malloc.c:236)
==2524==    by 0x408BF8: main (test1.c:291)
==2524== 
==2524== 1,296 bytes in 1 blocks are possibly lost in loss record 17 of 19
==2524==    at 0x4C274A8: malloc (vg_replace_malloc.c:236)
==2524==    by 0x404402: main (test1.c:93)
==2524== 
==2524== 1,296 bytes in 1 blocks are definitely lost in loss record 18 of 
19
==2524==    at 0x4C274A8: malloc (vg_replace_malloc.c:236)
==2524==    by 0x40437A: main (test1.c:74)
==2524== 
==2524== 1,296 bytes in 1 blocks are definitely lost in loss record 19 of 
19
==2524==    at 0x4C274A8: malloc (vg_replace_malloc.c:236)
==2524==    by 0x404387: main (test1.c:75)
==2524== 
==2524== LEAK SUMMARY:
==2524==    definitely lost: 3,456 bytes in 3 blocks
==2524==    indirectly lost: 0 bytes in 0 blocks
==2524==      possibly lost: 1,296 bytes in 1 blocks
==2524==    still reachable: 2,569 bytes in 15 blocks
==2524==         suppressed: 0 bytes in 0 blocks
==2524== 
==2524== For counts of detected and suppressed errors, rerun with: -v
==2524== ERROR SUMMARY: 5 errors from 5 contexts (suppressed: 4 from 4)
Segmentation fault
bruce at AEL-UBUS-01:~/sb/Private/Bruce_Labitt/ESPRIT_prod/C/esprit_proj/test_integration$ 


If the flag -fstack-protector-all is not set then valgrind output is the 
same. :(

Any thing I can try? 

-Bruce


******************************
Neither the footer nor anything else in this E-mail is intended to or constitutes an <br>electronic signature and/or legally binding agreement in the absence of an <br>express statement or Autoliv policy and/or procedure to the contrary.<br>This E-mail and any attachments hereto are Autoliv property and may contain legally <br>privileged, confidential and/or proprietary information.<br>The recipient of this E-mail is prohibited from distributing, copying, forwarding or in any way <br>disseminating any material contained within this E-mail without prior written <br>permission from the author. If you receive this E-mail in error, please <br>immediately notify the author and delete this E-mail.  Autoliv disclaims all <br>responsibility and liability for the consequences of any person who fails to <br>abide by the terms herein. <br>
******************************



More information about the gnhlug-discuss mailing list