LAPACK/ScaLAPACK Development

Hi all,
I have written dgesdd into some code and was using it to find the eigenvalues of quadratic matrices. For this example everything was ok. But if I use a matrix like 1397 x 106 entrys, Visual Studio crashes.

Here is a picture of the crash:


and here is my code:

Code: Select all

#include <time.h>
#include "stdio.h"
#include <iostream>

extern "C" {
   #include "f2c.h"
   #include "clapack.h"
}

using namespace std;

int main()
{
   char JOBU;
   char JOBVT;

   float startTime, endTime, actualTime;

   int i, j, k;

   int rows = 1397;
   int columns = 106;
   int size = rows * columns;

   integer M = rows;
   integer N = columns;

   integer LDA = M;
   integer LDU = M;
   integer LDVT = N;
   integer LWORK;
   integer IWORK;
   integer INFO;

   integer mn = min( M, N );
   integer MN = max( M, N );

   //double a[ROW*COLUMN] = {0,0,1,0,3,0,1,0,0};
   
   double *a = new double [size];
   for (int k = 0; k < size; k++)
      a[k] = 1;

   JOBU = 'A';
   JOBVT = 'A';

   LWORK = 56046;
   IWORK = 8*mn;

   double  *s      = new double[columns];
   double  *work   = new double[56046];
   integer *iwork  = new integer[IWORK];
   double  *uu     = new double[rows * columns];
   double  *vt     = new double[columns * columns];

   startTime = clock();

   //dgesvd_( &JOBU, &JOBVT, &M, &N, a, &LDA, s, uu, &LDU, vt, &LDVT, wk, &LWORK, &INFO);
   dgesdd_( &JOBU, &M, &N, a, &LDA, s, uu, &LDU, vt, &LDVT, work, &LWORK, iwork, &INFO);

   endTime = clock();
    actualTime = (endTime - startTime)/(CLOCKS_PER_SEC/1000);

   //printf("\n TIME=%f", actualTime );
   cout << "TIME: " << actualTime << "\n" << endl;
   
   //printf("\n INFO=%d", INFO );
   cout << "INFO: " << INFO << "\n" << endl;


   cout << "WORK: " << work[0] << "\n\n" << endl;

   cout << "s: " << s[0] << " " << s[1] << " " << s[2] << "\n\n" << endl;

   //for (i= 0; i< COLUMN; i++ )
   //{
   //   printf("\n s[ %d ] = %f", i, s[ i ] );
   //}

   //
   //for (j= 0; j< ROW*COLUMN; j++ )
   //{
   //      printf("\n uu[ %d ] = %f", j, uu[ j ] );
   //}

   //for (k= 0; k< ROW*COLUMN; k++ )
   //{
   //   printf("\n vt[ %d ] = %f", k, vt[ k ] );
   //}
   delete[] a;
   delete[] s;
   delete[] work;
   delete[] iwork;
   delete[] uu;
   delete[] vt;

   a     = NULL;
   s     = NULL;
   work  = NULL;
   iwork = NULL;
   uu    = NULL;
   vt    = NULL;

   return 0;
}

Does sombody have an idea what causes the error???

Hello

I look at your code.

I think the mistake is in JOBU. You can not set JOBU = 'A' and having U of size m-by-n.
JOBU='A' wil return the m-by-m unitary matrix U. JOBU = 'S' will return you only the first
n columns (case n<=m).

If you really want U of size m-by-n, then you should consider JOBU = 'O', This will save
you m*n memory space. The matrix A is in general useless after DGESDD anyhow so
why not storing U in A.

Below is a working code, as you can see I am not a big fan of C++ ... but I am sure you'll fix this.

Julien.

Code: Select all

#include "stdio.h"

extern "C" void dgesdd_ ( char *JOBZ, int *M, int *N, double *A, int *LDA,
double *S, double *U, int *LDU, double *VT, int *LDVT, double *WORK,
int *LWORK, int *IWORK, int *INFO );

int main()
{
   char JOBU;
   int k;
   int M = 1397;
   int N = 106;
   int LWORK;
   int INFO;
   double a[M*N];
   double  s[N];
   double  work[56046];
   int iwork[8*N];
   double  uu[M*N];
   double  vt[N*N];

   for ( k = 0; k < M*N; k++)
      a[k] = 1;

   LWORK = 56046;

   JOBU = 'S';
   dgesdd_( &JOBU, &M, &N, a, &M, s, uu, &M, vt, &N, work, &LWORK, iwork, &INFO);

   printf("\n INFO=%d", INFO );

   return 0;
}



Julien Langou wrote:Hello

I look at your code.

I think the mistake is in JOBU. You can not set JOBU = 'A' and having U of size m-by-n.
JOBU='A' wil return the m-by-m unitary matrix U. JOBU = 'S' will return you only the first
n columns (case n<=m).

If you really want U of size m-by-n, then you should consider JOBU = 'O', This will save
you m*n memory space. The matrix A is in general useless after DGESDD anyhow so
why not storing U in A.

Below is a working code, as you can see I am not a big fan of C++ ... but I am sure you'll fix this.

Julien.

Code: Select all

#include "stdio.h"

extern "C" void dgesdd_ ( char *JOBZ, int *M, int *N, double *A, int *LDA,
double *S, double *U, int *LDU, double *VT, int *LDVT, double *WORK,
int *LWORK, int *IWORK, int *INFO );

int main()
{
   char JOBU;
   int k;
   int M = 1397;
   int N = 106;
   int LWORK;
   int INFO;
   double a[M*N];
   double  s[N];
   double  work[56046];
   int iwork[8*N];
   double  uu[M*N];
   double  vt[N*N];

   for ( k = 0; k < M*N; k++)
      a[k] = 1;

   LWORK = 56046;

   JOBU = 'S';
   dgesdd_( &JOBU, &M, &N, a, &M, s, uu, &M, vt, &N, work, &LWORK, iwork, &INFO);

   printf("\n INFO=%d", INFO );

   return 0;
}



Hi. I'm new in Clapack. This is the more simple code I've found. Thanks Julien. The problem is, how I compile it? Thank you