[Trilinos-Users] Problems with Teuchos_LAPACK::GESVD

[Kemelli C. Estacio-Hiroms]

Hi everyone

I am not sure if this is the forum where I should post my doubts with
Teuchos_LAPACK::GESVD, but I will give it a try (and my apologies if I
am in the wrong place :) )

I need to encapsulate the function Teuchos_LAPACK::GESVD, for singular
value decomposition of a matrix in a complex code, and as a first
step, I wrote a simple code for it.

However, the results I am getting are not right and even after a good
amount of looking around, I can't figure out if the bug is mine or at
Trilinos/Teuchos/Lapack.

Would anybody, please, give me some insights? I am attacing my simple
code here, together with the Makefile I am using.
Thanks a lot!

Kemelli.

-- 
Kemelli C. Estacio-Hiroms
Postdoctoral Researcher
Institute for Computational Engineering and Sciences
The University of Texas at Austin
(512) 299-8403
kemelli at ices.utexas.edu
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Makefile_svd_kemelli_violeta
Type: application/octet-stream
Size: 621 bytes
Desc: not available
Url : https://software.sandia.gov/pipermail/trilinos-users/attachments/20121121/491fee09/attachment.obj 
-------------- next part --------------
/*
  For this example, we work with the four by three matrix A[i,j] =
  n*i+j+1.  The output of the program should be

A:
    1	2	3	
    4	5	6	
    7	8	9	
    10	11	12	

S:
    25.4624	1.29066	2.21303e-15	

U:
    -0.140877	0.824714	0.540521	
    -0.343946	0.426264	-0.654711	
    -0.547016	0.0278135	-0.312142	
    -0.750086	-0.370637	0.426331	
V:    
    -0.504533	-0.574516	-0.644498	
    -0.760776	-0.0571405	0.646495	
    -0.408248	0.816497	-0.408248	

This values have been compared to with Matlab values. 
But the results using the Lapack routine GESVD don't match.  
Something is wrong.
*/


#include <cmath>
#include <stdio.h>
#include <fstream>
#include <Teuchos_LAPACK.hpp>


int main(int nargs, char *args[])
{
  double **A, *S, **U, **VT;
  
  int m=4; // rows
  int n=3; // columns
  int minmn, maxmn;

  if (m>=n) { minmn = n; maxmn = m; } else { minmn = m; maxmn = n; }
    
  // Generate the matrix for my equation

  A = new double*[m];
  
  for (int i=0; i<m; i++) 
    A[i] = new double[n];
  
  for (int i=0; i<m; i++) {
    for (int j=0; j<n; j++){
      A[i][j] = n*i+j+1;      
    }
  }

  printf("\nA:\n");
  for (int i=0; i<m; i++) {
    for (int j=0; j<n; j++) 
      printf(" %3.2f \t",A[i][j]);
    printf("\n");
  } 

  // Allocate the storage space for the output arrays and call the solver
  S = new double[minmn];
  U = new double*[m]; for (int i=0; i<m; i++) U[i] = new double[minmn];
  VT = new double*[minmn]; for (int i=0; i<minmn; i++) VT[i] = new double[n];

  //SDV========================================
  char jobu, jobvt;
  int lda, ldu, ldvt, lwork, info;
  double  *work;
  
  jobu = 'S';
  jobvt = 'S';
  lda = m; // The leading dimension of the matrix a.
  ldu = m; // Left singular vector matrix
  ldvt = minmn; // Right singular vector matrix
  
  lwork = 15*maxmn; // Set up the work array, larger than needed.
  work = new double[lwork];
  
  double *rwork;
  int aux1= 5*minmn+7, aux2= 2*maxmn+2*minmn+1;
  int aux_dim;
  
  if (aux1>=aux2) { aux_dim = minmn*aux1; } else {aux_dim = minmn*aux2; }
  
  //LRWORK >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1) //From Lapack documentation
  rwork = new double[aux_dim];
  
  Teuchos::LAPACK<int, double> lapack;
  lapack.GESVD(jobu,jobvt,m,n,&A[0][0],lda,&S[0],&U[0][0],ldu,&VT[0][0],ldvt,&work[0],lwork,&rwork[0],&info);

  //========================================

  // Output the results
  printf("\nS:\n");
  for (int i=0; i<minmn; i++)  printf(" %3.2f \t",S[i]);

  printf("\n\nU:\n");
  for (int i=0; i<m; i++) {
    for (int j=0; j<m; j++) printf(" %3.2f \t",U[i][j]); 
  printf("\n");
  } 
  
  printf("\nVT:\n");
  for (int i=0; i<n; i++) {
    for (int j=0; j<n; j++) printf(" %3.2f \t",VT[i][j]);
  printf("\n");
  }  
  
  if (info == 0) printf("\nSuccess.\n\n");
  else printf("\nSomething went wrong.\n\n");
    
  return 0;
}