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_zgesdd_work( int matrix_order, char jobz, lapack_int m,
00038 lapack_int n, lapack_complex_double* a,
00039 lapack_int lda, double* s,
00040 lapack_complex_double* u, lapack_int ldu,
00041 lapack_complex_double* vt, lapack_int ldvt,
00042 lapack_complex_double* work, lapack_int lwork,
00043 double* rwork, lapack_int* iwork )
00044 {
00045 lapack_int info = 0;
00046 if( matrix_order == LAPACK_COL_MAJOR ) {
00047
00048 LAPACK_zgesdd( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
00049 &lwork, rwork, iwork, &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( jobz, 'a' ) ||
00055 LAPACKE_lsame( jobz, 's' ) ||
00056 ( LAPACKE_lsame( jobz, 'o' ) && m<n) ) ? m : 1;
00057 lapack_int ncols_u = ( LAPACKE_lsame( jobz, 'a' ) ||
00058 ( LAPACKE_lsame( jobz, 'o' ) && m<n) ) ? m :
00059 ( LAPACKE_lsame( jobz, 's' ) ? MIN(m,n) : 1);
00060 lapack_int nrows_vt = ( LAPACKE_lsame( jobz, 'a' ) ||
00061 ( LAPACKE_lsame( jobz, 'o' ) && m<n) ) ? n :
00062 ( LAPACKE_lsame( jobz, 's' ) ? MIN(m,n) : 1);
00063 lapack_int lda_t = MAX(1,m);
00064 lapack_int ldu_t = MAX(1,nrows_u);
00065 lapack_int ldvt_t = MAX(1,nrows_vt);
00066 lapack_complex_double* a_t = NULL;
00067 lapack_complex_double* u_t = NULL;
00068 lapack_complex_double* vt_t = NULL;
00069
00070 if( lda < n ) {
00071 info = -6;
00072 LAPACKE_xerbla( "LAPACKE_zgesdd_work", info );
00073 return info;
00074 }
00075 if( ldu < ncols_u ) {
00076 info = -9;
00077 LAPACKE_xerbla( "LAPACKE_zgesdd_work", info );
00078 return info;
00079 }
00080 if( ldvt < n ) {
00081 info = -11;
00082 LAPACKE_xerbla( "LAPACKE_zgesdd_work", info );
00083 return info;
00084 }
00085
00086 if( lwork == -1 ) {
00087 LAPACK_zgesdd( &jobz, &m, &n, a, &lda_t, s, u, &ldu_t, vt, &ldvt_t,
00088 work, &lwork, rwork, iwork, &info );
00089 return (info < 0) ? (info - 1) : info;
00090 }
00091
00092 a_t = (lapack_complex_double*)
00093 LAPACKE_malloc( sizeof(lapack_complex_double) * lda_t * MAX(1,n) );
00094 if( a_t == NULL ) {
00095 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00096 goto exit_level_0;
00097 }
00098 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00099 ( LAPACKE_lsame( jobz, 'o' ) && (m<n) ) ) {
00100 u_t = (lapack_complex_double*)
00101 LAPACKE_malloc( sizeof(lapack_complex_double) *
00102 ldu_t * MAX(1,ncols_u) );
00103 if( u_t == NULL ) {
00104 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00105 goto exit_level_1;
00106 }
00107 }
00108 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00109 ( LAPACKE_lsame( jobz, 'o' ) && (m>=n) ) ) {
00110 vt_t = (lapack_complex_double*)
00111 LAPACKE_malloc( sizeof(lapack_complex_double) *
00112 ldvt_t * MAX(1,n) );
00113 if( vt_t == NULL ) {
00114 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00115 goto exit_level_2;
00116 }
00117 }
00118
00119 LAPACKE_zge_trans( matrix_order, m, n, a, lda, a_t, lda_t );
00120
00121 LAPACK_zgesdd( &jobz, &m, &n, a_t, &lda_t, s, u_t, &ldu_t, vt_t,
00122 &ldvt_t, work, &lwork, rwork, iwork, &info );
00123 if( info < 0 ) {
00124 info = info - 1;
00125 }
00126
00127 LAPACKE_zge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
00128 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00129 ( LAPACKE_lsame( jobz, 'o' ) && (m<n) ) ) {
00130 LAPACKE_zge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
00131 u, ldu );
00132 }
00133 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00134 ( LAPACKE_lsame( jobz, 'o' ) && (m>=n) ) ) {
00135 LAPACKE_zge_trans( LAPACK_COL_MAJOR, nrows_vt, n, vt_t, ldvt_t, vt,
00136 ldvt );
00137 }
00138
00139 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00140 ( LAPACKE_lsame( jobz, 'o' ) && (m>=n) ) ) {
00141 LAPACKE_free( vt_t );
00142 }
00143 exit_level_2:
00144 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00145 ( LAPACKE_lsame( jobz, 'o' ) && (m<n) ) ) {
00146 LAPACKE_free( u_t );
00147 }
00148 exit_level_1:
00149 LAPACKE_free( a_t );
00150 exit_level_0:
00151 if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
00152 LAPACKE_xerbla( "LAPACKE_zgesdd_work", info );
00153 }
00154 } else {
00155 info = -1;
00156 LAPACKE_xerbla( "LAPACKE_zgesdd_work", info );
00157 }
00158 return info;
00159 }