/*  svd.c  (C) 1997 Mark A. Livingston, Shankar Krishnan */
/* Any non-commercial use of this code is permitted so long as */
/* the above copyright line appears in the file.  Please contact */
/* the author regarding permission for commercial use. */
#include <stdio.h>
#include <stdlib.h>
#include <quat.h>
#include "svd.h"

/* compute the singular value decomposition of a real m x n matrix mat */
/* this assumes u2, sigma, v2 have been alloc'ed before calling this routine */
/* mat is unchanged at the end of the call */
void
call_svd( q_matrix_type mat, int m, int n, q_matrix_type rgrgu,
	  q_vec_type rgsigma, q_matrix_type rgrgv )
{ 
    int     i, j, k, ierr;
    int     min_mn, lda, ldu, lwork, ldvt;
    double *a, *u, *vt, *s, *e, *work;
    char    jobu, jobvt;
    int     job;
    double  LAmag, LAsum;
    int     LAi, LAj, LAk;

    /* allocate/initialize the matrices */
    a = (double *) malloc( ( m * n ) * sizeof( double ) );
    for( j = k = 0; j < n; j++ )
      for( i = 0; i < m; i++ )
	a[k++] = mat[i][j];	

    min_mn = ( m > n ) ? n : m;

    /*
       jobu = 'A';
       jobvt = 'A';
    */
    job = 11;
    lda = ldu = m;
    ldvt = n;

    lwork = ( m > n ) ? 5 * m : 5 * n;
    u    = (double *) malloc( ( m * m ) * sizeof( double ) );
    vt   = (double *) malloc( ( n * n ) * sizeof( double ) );
    work = (double *) malloc( ( lwork ) * sizeof( double ) );
    s    = (double *) malloc( ( min_mn ) * sizeof( double ) );
    e    = (double *) malloc( ( n ) * sizeof( double ) );

    /* call the svd subroutine */
    /*
    dgesvd_( &jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, 
             &lwork, &ierr );
    */
#ifdef __SGI
    dsvdc_( a, &lda, &m, &n, s, e, u, &ldu, vt, &ldvt, work, &job, &ierr );
#else
    dsvdc( a, &lda, &m, &n, s, e, u, &ldu, vt, &ldvt, work, &job, &ierr );
#endif

    if( ierr != 0 ) {
      printf( "Error computing SVD\n" ); 
    }

    /* form the 2-dim arrays u, v, and the 1-dim array sigma */
    for( k = j = 0; j < m; j++ )
      for( i = 0; i < m; i++ )
	rgrgu[i][j] = u[k++];

    for( k = j = 0; j < n; j++ )
      for( i = 0; i < n; i++ )
	rgrgv[i][j] = vt[k++]; /* This was changed from [j][i] to [i][j]
				  when I changed from dgesvd to dsvdc */

    for( LAi = 0; LAi < min_mn; LAi++ )
      rgsigma[ LAi ] = s[ LAi ];
/*     VEC_ASN_OP( rgsigma, =, s, min_mn ); */

    /* free up the allocated space */
    free( a );
    free( u );
    free( vt );
    free( work );
    free( s );

    return;
}

void
print_fort_mat( double *mat, int m, int n )
{ 
    int	i, j, k; 

    for( i = 0; i < m; i++ ) {
      k = i; 
      for( j = 0; j < n; j++ ) {
	printf( "%f ", mat[k] );
	k += m;
      }
      printf( "\n" );
    }
    return;
}