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