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_sgesdd_work( int matrix_order, char jobz, lapack_int m,
00038 lapack_int n, float* a, lapack_int lda,
00039 float* s, float* u, lapack_int ldu, float* vt,
00040 lapack_int ldvt, float* work, lapack_int lwork,
00041 lapack_int* iwork )
00042 {
00043 lapack_int info = 0;
00044 if( matrix_order == LAPACK_COL_MAJOR ) {
00045
00046 LAPACK_sgesdd( &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 float* a_t = NULL;
00065 float* u_t = NULL;
00066 float* vt_t = NULL;
00067
00068 if( lda < n ) {
00069 info = -6;
00070 LAPACKE_xerbla( "LAPACKE_sgesdd_work", info );
00071 return info;
00072 }
00073 if( ldu < ncols_u ) {
00074 info = -9;
00075 LAPACKE_xerbla( "LAPACKE_sgesdd_work", info );
00076 return info;
00077 }
00078 if( ldvt < n ) {
00079 info = -11;
00080 LAPACKE_xerbla( "LAPACKE_sgesdd_work", info );
00081 return info;
00082 }
00083
00084 if( lwork == -1 ) {
00085 LAPACK_sgesdd( &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 = (float*)LAPACKE_malloc( sizeof(float) * 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 = (float*)
00098 LAPACKE_malloc( sizeof(float) * 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 = (float*)LAPACKE_malloc( sizeof(float) * ldvt_t * MAX(1,n) );
00107 if( vt_t == NULL ) {
00108 info = LAPACK_TRANSPOSE_MEMORY_ERROR;
00109 goto exit_level_2;
00110 }
00111 }
00112
00113 LAPACKE_sge_trans( matrix_order, m, n, a, lda, a_t, lda_t );
00114
00115 LAPACK_sgesdd( &jobz, &m, &n, a_t, &lda_t, s, u_t, &ldu_t, vt_t,
00116 &ldvt_t, work, &lwork, iwork, &info );
00117 if( info < 0 ) {
00118 info = info - 1;
00119 }
00120
00121 LAPACKE_sge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
00122 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00123 ( LAPACKE_lsame( jobz, 'o' ) && (m<n) ) ) {
00124 LAPACKE_sge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
00125 u, ldu );
00126 }
00127 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00128 ( LAPACKE_lsame( jobz, 'o' ) && (m>=n) ) ) {
00129 LAPACKE_sge_trans( LAPACK_COL_MAJOR, nrows_vt, n, vt_t, ldvt_t, vt,
00130 ldvt );
00131 }
00132
00133 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00134 ( LAPACKE_lsame( jobz, 'o' ) && (m>=n) ) ) {
00135 LAPACKE_free( vt_t );
00136 }
00137 exit_level_2:
00138 if( LAPACKE_lsame( jobz, 'a' ) || LAPACKE_lsame( jobz, 's' ) ||
00139 ( LAPACKE_lsame( jobz, 'o' ) && (m<n) ) ) {
00140 LAPACKE_free( u_t );
00141 }
00142 exit_level_1:
00143 LAPACKE_free( a_t );
00144 exit_level_0:
00145 if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
00146 LAPACKE_xerbla( "LAPACKE_sgesdd_work", info );
00147 }
00148 } else {
00149 info = -1;
00150 LAPACKE_xerbla( "LAPACKE_sgesdd_work", info );
00151 }
00152 return info;
00153 }