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_sgesvd_work( int matrix_order, char jobu, char jobvt,
00038 lapack_int m, lapack_int n, float* a,
00039 lapack_int lda, float* s, float* u,
00040 lapack_int ldu, float* vt, lapack_int ldvt,
00041 float* work, lapack_int lwork )
00042 {
00043 lapack_int info = 0;
00044 if( matrix_order == LAPACK_COL_MAJOR ) {
00045
00046 LAPACK_sgesvd( &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 float* a_t = NULL;
00062 float* u_t = NULL;
00063 float* vt_t = NULL;
00064
00065 if( lda < n ) {
00066 info = -7;
00067 LAPACKE_xerbla( "LAPACKE_sgesvd_work", info );
00068 return info;
00069 }
00070 if( ldu < ncols_u ) {
00071 info = -10;
00072 LAPACKE_xerbla( "LAPACKE_sgesvd_work", info );
00073 return info;
00074 }
00075 if( ldvt < n ) {
00076 info = -12;
00077 LAPACKE_xerbla( "LAPACKE_sgesvd_work", info );
00078 return info;
00079 }
00080
00081 if( lwork == -1 ) {
00082 LAPACK_sgesvd( &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 = (float*)LAPACKE_malloc( sizeof(float) * 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 = (float*)
00094 LAPACKE_malloc( sizeof(float) * 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 = (float*)LAPACKE_malloc( sizeof(float) * ldvt_t * MAX(1,n) );
00102 if( vt_t == NULL ) {
00103 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00104 goto exit_level_2;
00105 }
00106 }
00107
00108 LAPACKE_sge_trans( matrix_order, m, n, a, lda, a_t, lda_t );
00109
00110 LAPACK_sgesvd( &jobu, &jobvt, &m, &n, a_t, &lda_t, s, u_t, &ldu_t, vt_t,
00111 &ldvt_t, work, &lwork, &info );
00112 if( info < 0 ) {
00113 info = info - 1;
00114 }
00115
00116 LAPACKE_sge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
00117 if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
00118 LAPACKE_sge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
00119 u, ldu );
00120 }
00121 if( LAPACKE_lsame( jobvt, 'a' ) || LAPACKE_lsame( jobvt, 's' ) ) {
00122 LAPACKE_sge_trans( LAPACK_COL_MAJOR, nrows_vt, n, vt_t, ldvt_t, vt,
00123 ldvt );
00124 }
00125
00126 if( LAPACKE_lsame( jobvt, 'a' ) || LAPACKE_lsame( jobvt, 's' ) ) {
00127 LAPACKE_free( vt_t );
00128 }
00129 exit_level_2:
00130 if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
00131 LAPACKE_free( u_t );
00132 }
00133 exit_level_1:
00134 LAPACKE_free( a_t );
00135 exit_level_0:
00136 if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
00137 LAPACKE_xerbla( "LAPACKE_sgesvd_work", info );
00138 }
00139 } else {
00140 info = -1;
00141 LAPACKE_xerbla( "LAPACKE_sgesvd_work", info );
00142 }
00143 return info;
00144 }