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