libflame  revision_anchor
Functions
FLA_Svd_external.c File Reference

(r)

Functions

FLA_Error FLA_Svd_external (FLA_Svd_type jobu, FLA_Svd_type jobv, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)
 

Function Documentation

◆ FLA_Svd_external()

FLA_Error FLA_Svd_external ( FLA_Svd_type  jobu,
FLA_Svd_type  jobv,
FLA_Obj  A,
FLA_Obj  s,
FLA_Obj  U,
FLA_Obj  V 
)

References F77_cgesvd(), F77_dgesvd(), F77_sgesvd(), F77_zgesvd(), FLA_Check_error_level(), FLA_Obj_col_stride(), FLA_Obj_create(), FLA_Obj_datatype(), FLA_Obj_datatype_proj_to_real(), FLA_Obj_free(), FLA_Obj_has_zero_dim(), FLA_Obj_is_complex(), FLA_Obj_length(), FLA_Obj_width(), FLA_Param_map_flame_to_netlib_svd_type(), FLA_Svd_check(), and i.

14 {
15  int info = 0;
16 #ifdef FLA_ENABLE_EXTERNAL_LAPACK_INTERFACES
17  FLA_Datatype datatype;
18  FLA_Datatype dt_real;
19  int m_A, n_A, cs_A;
20  int cs_U;
21  int cs_V;
22  int min_m_n;
23  int lwork, lrwork;
24  FLA_Obj work, rwork;
25  char blas_jobu;
26  char blas_jobv;
27  int i;
28 
29  if ( FLA_Check_error_level() == FLA_FULL_ERROR_CHECKING )
30  FLA_Svd_check( jobu, jobv, A, s, U, V );
31 
32  if ( FLA_Obj_has_zero_dim( A ) ) return FLA_SUCCESS;
33 
34  datatype = FLA_Obj_datatype( A );
35  dt_real = FLA_Obj_datatype_proj_to_real( A );
36 
37  m_A = FLA_Obj_length( A );
38  n_A = FLA_Obj_width( A );
39  cs_A = FLA_Obj_col_stride( A );
40 
41  cs_U = FLA_Obj_col_stride( U );
42 
43  cs_V = FLA_Obj_col_stride( V );
44 
45  min_m_n = min( m_A, n_A );
46 
47  // Allocate the rwork array up front since its size is not dependent on
48  // internal block sizes.
49  lrwork = 5 * min_m_n;
50  if ( FLA_Obj_is_complex( A ) )
51  FLA_Obj_create( dt_real, lrwork, 1, 0, 0, &rwork );
52 
53  FLA_Param_map_flame_to_netlib_svd_type( jobu, &blas_jobu );
54  FLA_Param_map_flame_to_netlib_svd_type( jobv, &blas_jobv );
55 
56  // Make a workspace query the first time through. This will provide us with
57  // and ideal workspace size based on an internal block size.
58  lwork = -1;
59  FLA_Obj_create( datatype, 1, 1, 0, 0, &work );
60 
61  for ( i = 0; i < 2; ++i )
62  {
63  if ( i == 1 )
64  {
65  // Grab the queried ideal workspace size from the work array, free the
66  // work object, and then re-allocate the workspace with the ideal size.
67  if ( datatype == FLA_FLOAT || datatype == FLA_COMPLEX )
68  lwork = ( int ) *FLA_FLOAT_PTR( work );
69  else if ( datatype == FLA_DOUBLE || datatype == FLA_DOUBLE_COMPLEX )
70  lwork = ( int ) *FLA_DOUBLE_PTR( work );
71 
72  FLA_Obj_free( &work );
73  FLA_Obj_create( datatype, lwork, 1, 0, 0, &work );
74  }
75 
76  switch( datatype ) {
77 
78  case FLA_FLOAT:
79  {
80  float* buff_A = ( float* ) FLA_FLOAT_PTR( A );
81  float* buff_s = ( float* ) FLA_FLOAT_PTR( s );
82  float* buff_U = ( float* ) FLA_FLOAT_PTR( U );
83  float* buff_V = ( float* ) FLA_FLOAT_PTR( V );
84  float* buff_work = ( float* ) FLA_FLOAT_PTR( work );
85 
86  F77_sgesvd( &blas_jobu,
87  &blas_jobv,
88  &m_A,
89  &n_A,
90  buff_A, &cs_A,
91  buff_s,
92  buff_U, &cs_U,
93  buff_V, &cs_V,
94  buff_work, &lwork,
95  &info );
96 
97  break;
98  }
99 
100  case FLA_DOUBLE:
101  {
102  double* buff_A = ( double* ) FLA_DOUBLE_PTR( A );
103  double* buff_s = ( double* ) FLA_DOUBLE_PTR( s );
104  double* buff_U = ( double* ) FLA_DOUBLE_PTR( U );
105  double* buff_V = ( double* ) FLA_DOUBLE_PTR( V );
106  double* buff_work = ( double* ) FLA_DOUBLE_PTR( work );
107 
108  F77_dgesvd( &blas_jobu,
109  &blas_jobv,
110  &m_A,
111  &n_A,
112  buff_A, &cs_A,
113  buff_s,
114  buff_U, &cs_U,
115  buff_V, &cs_V,
116  buff_work, &lwork,
117  &info );
118 
119  break;
120  }
121 
122  case FLA_COMPLEX:
123  {
124  scomplex* buff_A = ( scomplex* ) FLA_COMPLEX_PTR( A );
125  float* buff_s = ( float* ) FLA_FLOAT_PTR( s );
126  scomplex* buff_U = ( scomplex* ) FLA_COMPLEX_PTR( U );
127  scomplex* buff_V = ( scomplex* ) FLA_COMPLEX_PTR( V );
128  scomplex* buff_work = ( scomplex* ) FLA_COMPLEX_PTR( work );
129  float* buff_rwork = ( float* ) FLA_FLOAT_PTR( rwork );
130 
131  F77_cgesvd( &blas_jobu,
132  &blas_jobv,
133  &m_A,
134  &n_A,
135  buff_A, &cs_A,
136  buff_s,
137  buff_U, &cs_U,
138  buff_V, &cs_V,
139  buff_work, &lwork,
140  buff_rwork,
141  &info );
142 
143  break;
144  }
145 
146  case FLA_DOUBLE_COMPLEX:
147  {
148  dcomplex* buff_A = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( A );
149  double* buff_s = ( double* ) FLA_DOUBLE_PTR( s );
150  dcomplex* buff_U = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( U );
151  dcomplex* buff_V = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( V );
152  dcomplex* buff_work = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( work );
153  double* buff_rwork = ( double* ) FLA_DOUBLE_PTR( rwork );
154 
155  F77_zgesvd( &blas_jobu,
156  &blas_jobv,
157  &m_A,
158  &n_A,
159  buff_A, &cs_A,
160  buff_s,
161  buff_U, &cs_U,
162  buff_V, &cs_V,
163  buff_work, &lwork,
164  buff_rwork,
165  &info );
166 
167  break;
168  }
169 
170  }
171  }
172 
173  FLA_Obj_free( &work );
174  if ( FLA_Obj_is_complex( A ) )
175  FLA_Obj_free( &rwork );
176 #else
177  FLA_Check_error_code( FLA_EXTERNAL_LAPACK_NOT_IMPLEMENTED );
178 #endif
179 
180  return info;
181 }
FLA_Error FLA_Obj_create(FLA_Datatype datatype, dim_t m, dim_t n, dim_t rs, dim_t cs, FLA_Obj *obj)
Definition: FLA_Obj.c:55
FLA_Error FLA_Obj_free(FLA_Obj *obj)
Definition: FLA_Obj.c:588
int F77_zgesvd(char *jobu, char *jobv, int *m, int *n, dcomplex *a, int *lda, double *s, dcomplex *u, int *ldu, dcomplex *vt, int *ldvt, dcomplex *work, int *lwork, double *rwork, int *info)
FLA_Datatype FLA_Obj_datatype(FLA_Obj obj)
Definition: FLA_Query.c:13
Definition: FLA_type_defs.h:158
FLA_Bool FLA_Obj_has_zero_dim(FLA_Obj A)
Definition: FLA_Query.c:400
dim_t FLA_Obj_width(FLA_Obj obj)
Definition: FLA_Query.c:123
int F77_sgesvd(char *jobu, char *jobv, int *m, int *n, float *a, int *lda, float *s, float *u, int *ldu, float *vt, int *ldvt, float *work, int *lwork, int *info)
FLA_Datatype FLA_Obj_datatype_proj_to_real(FLA_Obj A)
Definition: FLA_Query.c:23
Definition: blis_type_defs.h:132
void FLA_Param_map_flame_to_netlib_svd_type(FLA_Svd_type svd_type, void *lapack_svd_type)
Definition: FLA_Param.c:171
FLA_Error FLA_Svd_check(FLA_Svd_type jobu, FLA_Svd_type jobv, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)
Definition: FLA_Svd_check.c:13
unsigned int FLA_Check_error_level(void)
Definition: FLA_Check.c:18
int FLA_Datatype
Definition: FLA_type_defs.h:49
dim_t FLA_Obj_col_stride(FLA_Obj obj)
Definition: FLA_Query.c:174
int i
Definition: bl1_axmyv2.c:145
FLA_Bool FLA_Obj_is_complex(FLA_Obj A)
Definition: FLA_Query.c:324
dim_t FLA_Obj_length(FLA_Obj obj)
Definition: FLA_Query.c:116
int F77_dgesvd(char *jobu, char *jobv, int *m, int *n, double *a, int *lda, double *s, double *u, int *ldu, double *vt, int *ldvt, double *work, int *lwork, int *info)
Definition: blis_type_defs.h:137
int F77_cgesvd(char *jobu, char *jobv, int *m, int *n, scomplex *a, int *lda, float *s, scomplex *u, int *ldu, scomplex *vt, int *ldvt, scomplex *work, int *lwork, float *rwork, int *info)