00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "lapacke.h"
00035 #include "lapacke_utils.h"
00036
00037 lapack_int LAPACKE_zgesvd_work( int matrix_order, char jobu, char jobvt,
00038 lapack_int m, lapack_int n,
00039 lapack_complex_double* a, lapack_int lda,
00040 double* s, lapack_complex_double* u,
00041 lapack_int ldu, lapack_complex_double* vt,
00042 lapack_int ldvt, lapack_complex_double* work,
00043 lapack_int lwork, double* rwork )
00044 {
00045 lapack_int info = 0;
00046 if( matrix_order == LAPACK_COL_MAJOR ) {
00047
00048 LAPACK_zgesvd( &jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
00049 work, &lwork, rwork, &info );
00050 if( info < 0 ) {
00051 info = info - 1;
00052 }
00053 } else if( matrix_order == LAPACK_ROW_MAJOR ) {
00054 lapack_int nrows_u = ( LAPACKE_lsame( jobu, 'a' ) ||
00055 LAPACKE_lsame( jobu, 's' ) ) ? m : 1;
00056 lapack_int ncols_u = LAPACKE_lsame( jobu, 'a' ) ? m :
00057 ( LAPACKE_lsame( jobu, 's' ) ? MIN(m,n) : 1);
00058 lapack_int nrows_vt = LAPACKE_lsame( jobvt, 'a' ) ? n :
00059 ( LAPACKE_lsame( jobvt, 's' ) ? MIN(m,n) : 1);
00060 lapack_int lda_t = MAX(1,m);
00061 lapack_int ldu_t = MAX(1,nrows_u);
00062 lapack_int ldvt_t = MAX(1,nrows_vt);
00063 lapack_complex_double* a_t = NULL;
00064 lapack_complex_double* u_t = NULL;
00065 lapack_complex_double* vt_t = NULL;
00066
00067 if( lda < n ) {
00068 info = -7;
00069 LAPACKE_xerbla( "LAPACKE_zgesvd_work", info );
00070 return info;
00071 }
00072 if( ldu < ncols_u ) {
00073 info = -10;
00074 LAPACKE_xerbla( "LAPACKE_zgesvd_work", info );
00075 return info;
00076 }
00077 if( ldvt < n ) {
00078 info = -12;
00079 LAPACKE_xerbla( "LAPACKE_zgesvd_work", info );
00080 return info;
00081 }
00082
00083 if( lwork == -1 ) {
00084 LAPACK_zgesvd( &jobu, &jobvt, &m, &n, a, &lda_t, s, u, &ldu_t, vt,
00085 &ldvt_t, work, &lwork, rwork, &info );
00086 return (info < 0) ? (info - 1) : info;
00087 }
00088
00089 a_t = (lapack_complex_double*)
00090 LAPACKE_malloc( sizeof(lapack_complex_double) * lda_t * MAX(1,n) );
00091 if( a_t == NULL ) {
00092 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00093 goto exit_level_0;
00094 }
00095 if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
00096 u_t = (lapack_complex_double*)
00097 LAPACKE_malloc( sizeof(lapack_complex_double) *
00098 ldu_t * MAX(1,ncols_u) );
00099 if( u_t == NULL ) {
00100 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00101 goto exit_level_1;
00102 }
00103 }
00104 if( LAPACKE_lsame( jobvt, 'a' ) || LAPACKE_lsame( jobvt, 's' ) ) {
00105 vt_t = (lapack_complex_double*)
00106 LAPACKE_malloc( sizeof(lapack_complex_double) *
00107 ldvt_t * MAX(1,n) );
00108 if( vt_t == NULL ) {
00109 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00110 goto exit_level_2;
00111 }
00112 }
00113
00114 LAPACKE_zge_trans( matrix_order, m, n, a, lda, a_t, lda_t );
00115
00116 LAPACK_zgesvd( &jobu, &jobvt, &m, &n, a_t, &lda_t, s, u_t, &ldu_t, vt_t,
00117 &ldvt_t, work, &lwork, rwork, &info );
00118 if( info < 0 ) {
00119 info = info - 1;
00120 }
00121
00122 LAPACKE_zge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
00123 if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
00124 LAPACKE_zge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
00125 u, ldu );
00126 }
00127 if( LAPACKE_lsame( jobvt, 'a' ) || LAPACKE_lsame( jobvt, 's' ) ) {
00128 LAPACKE_zge_trans( LAPACK_COL_MAJOR, nrows_vt, n, vt_t, ldvt_t, vt,
00129 ldvt );
00130 }
00131
00132 if( LAPACKE_lsame( jobvt, 'a' ) || LAPACKE_lsame( jobvt, 's' ) ) {
00133 LAPACKE_free( vt_t );
00134 }
00135 exit_level_2:
00136 if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
00137 LAPACKE_free( u_t );
00138 }
00139 exit_level_1:
00140 LAPACKE_free( a_t );
00141 exit_level_0:
00142 if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
00143 LAPACKE_xerbla( "LAPACKE_zgesvd_work", info );
00144 }
00145 } else {
00146 info = -1;
00147 LAPACKE_xerbla( "LAPACKE_zgesvd_work", info );
00148 }
00149 return info;
00150 }