libflame  revision_anchor
Functions
FLA_Svdd_external.c File Reference

(r)

Functions

FLA_Error FLA_Svdd_external (FLA_Svd_type jobz, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)
 

Function Documentation

◆ FLA_Svdd_external()

FLA_Error FLA_Svdd_external ( FLA_Svd_type  jobz,
FLA_Obj  A,
FLA_Obj  s,
FLA_Obj  U,
FLA_Obj  V 
)

References F77_cgesdd(), F77_dgesdd(), F77_sgesdd(), F77_zgesdd(), 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_Svdd_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  FLA_Datatype dt_int;
20  int m_A, n_A, cs_A;
21  int cs_U;
22  int cs_V;
23  int min_m_n;
24  int lwork, lrwork, liwork;
25  FLA_Obj work, rwork, iwork;
26  char blas_jobz;
27  int i;
28 
29  if ( FLA_Check_error_level() == FLA_FULL_ERROR_CHECKING )
30  FLA_Svdd_check( jobz, 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  dt_int = FLA_INT;
37 
38  m_A = FLA_Obj_length( A );
39  n_A = FLA_Obj_width( A );
40  cs_A = FLA_Obj_col_stride( A );
41 
42  cs_U = FLA_Obj_col_stride( U );
43 
44  cs_V = FLA_Obj_col_stride( V );
45 
46  min_m_n = min( m_A, n_A );
47 
48  // Allocate the rwork and iwork arrays up front.
49  if ( jobz == FLA_SVD_VECTORS_NONE ) lrwork = 5 * min_m_n;
50  else lrwork = 5 * min_m_n * min_m_n +
51  7 * min_m_n;
52  liwork = 8 * min_m_n;
53 
54  FLA_Obj_create( dt_int, liwork, 1, 0, 0, &iwork );
55  if ( FLA_Obj_is_complex( A ) )
56  FLA_Obj_create( dt_real, lrwork, 1, 0, 0, &rwork );
57 
58  FLA_Param_map_flame_to_netlib_svd_type( jobz, &blas_jobz );
59 
60  // Make a workspace query the first time through. This will provide us with
61  // and ideal workspace size based on an internal block size.
62  lwork = -1;
63  FLA_Obj_create( datatype, 1, 1, 0, 0, &work );
64 
65  for ( i = 0; i < 2; ++i )
66  {
67  if ( i == 1 )
68  {
69  // Grab the queried ideal workspace size from the work array, free the
70  // work object, and then re-allocate the workspace with the ideal size.
71  if ( datatype == FLA_FLOAT || datatype == FLA_COMPLEX )
72  lwork = ( int ) *FLA_FLOAT_PTR( work );
73  else if ( datatype == FLA_DOUBLE || datatype == FLA_DOUBLE_COMPLEX )
74  lwork = ( int ) *FLA_DOUBLE_PTR( work );
75 
76  FLA_Obj_free( &work );
77  FLA_Obj_create( datatype, lwork, 1, 0, 0, &work );
78  }
79 
80  switch( datatype ) {
81 
82  case FLA_FLOAT:
83  {
84  float* buff_A = ( float* ) FLA_FLOAT_PTR( A );
85  float* buff_s = ( float* ) FLA_FLOAT_PTR( s );
86  float* buff_U = ( float* ) FLA_FLOAT_PTR( U );
87  float* buff_V = ( float* ) FLA_FLOAT_PTR( V );
88  float* buff_work = ( float* ) FLA_FLOAT_PTR( work );
89  int* buff_iwork = ( int* ) FLA_INT_PTR( iwork );
90 
91  F77_sgesdd( &blas_jobz,
92  &m_A,
93  &n_A,
94  buff_A, &cs_A,
95  buff_s,
96  buff_U, &cs_U,
97  buff_V, &cs_V,
98  buff_work, &lwork,
99  buff_iwork,
100  &info );
101 
102  break;
103  }
104 
105  case FLA_DOUBLE:
106  {
107  double* buff_A = ( double* ) FLA_DOUBLE_PTR( A );
108  double* buff_s = ( double* ) FLA_DOUBLE_PTR( s );
109  double* buff_U = ( double* ) FLA_DOUBLE_PTR( U );
110  double* buff_V = ( double* ) FLA_DOUBLE_PTR( V );
111  double* buff_work = ( double* ) FLA_DOUBLE_PTR( work );
112  int* buff_iwork = ( int* ) FLA_INT_PTR( iwork );
113 
114  F77_dgesdd( &blas_jobz,
115  &m_A,
116  &n_A,
117  buff_A, &cs_A,
118  buff_s,
119  buff_U, &cs_U,
120  buff_V, &cs_V,
121  buff_work, &lwork,
122  buff_iwork,
123  &info );
124 
125  break;
126  }
127 
128  case FLA_COMPLEX:
129  {
130  scomplex* buff_A = ( scomplex* ) FLA_COMPLEX_PTR( A );
131  float* buff_s = ( float* ) FLA_FLOAT_PTR( s );
132  scomplex* buff_U = ( scomplex* ) FLA_COMPLEX_PTR( U );
133  scomplex* buff_V = ( scomplex* ) FLA_COMPLEX_PTR( V );
134  scomplex* buff_work = ( scomplex* ) FLA_COMPLEX_PTR( work );
135  float* buff_rwork = ( float* ) FLA_FLOAT_PTR( rwork );
136  int* buff_iwork = ( int* ) FLA_INT_PTR( iwork );
137 
138  F77_cgesdd( &blas_jobz,
139  &m_A,
140  &n_A,
141  buff_A, &cs_A,
142  buff_s,
143  buff_U, &cs_U,
144  buff_V, &cs_V,
145  buff_work, &lwork,
146  buff_rwork,
147  buff_iwork,
148  &info );
149 
150  break;
151  }
152 
153  case FLA_DOUBLE_COMPLEX:
154  {
155  dcomplex* buff_A = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( A );
156  double* buff_s = ( double* ) FLA_DOUBLE_PTR( s );
157  dcomplex* buff_U = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( U );
158  dcomplex* buff_V = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( V );
159  dcomplex* buff_work = ( dcomplex* ) FLA_DOUBLE_COMPLEX_PTR( work );
160  double* buff_rwork = ( double* ) FLA_DOUBLE_PTR( rwork );
161  int* buff_iwork = ( int* ) FLA_INT_PTR( iwork );
162 
163  F77_zgesdd( &blas_jobz,
164  &m_A,
165  &n_A,
166  buff_A, &cs_A,
167  buff_s,
168  buff_U, &cs_U,
169  buff_V, &cs_V,
170  buff_work, &lwork,
171  buff_rwork,
172  buff_iwork,
173  &info );
174 
175  break;
176  }
177 
178  }
179  }
180 
181  FLA_Obj_free( &work );
182  FLA_Obj_free( &iwork );
183  if ( FLA_Obj_is_complex( A ) )
184  FLA_Obj_free( &rwork );
185 #else
186  FLA_Check_error_code( FLA_EXTERNAL_LAPACK_NOT_IMPLEMENTED );
187 #endif
188 
189  return info;
190 }
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
int F77_dgesdd(char *jobz, int *m, int *n, double *a, int *lda, double *s, double *u, int *ldu, double *vt, int *ldvt, double *work, int *lwork, int *iwork, int *info)
FLA_Error FLA_Obj_free(FLA_Obj *obj)
Definition: FLA_Obj.c:588
int F77_cgesdd(char *jobz, 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 *iwork, int *info)
int F77_zgesdd(char *jobz, 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 *iwork, 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_sgesdd(char *jobz, int *m, int *n, float *a, int *lda, float *s, float *u, int *ldu, float *vt, int *ldvt, float *work, int *lwork, int *iwork, 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
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
FLA_Error FLA_Svdd_check(FLA_Svd_type jobz, FLA_Obj A, FLA_Obj s, FLA_Obj U, FLA_Obj V)
Definition: FLA_Svdd_check.c:13
dim_t FLA_Obj_length(FLA_Obj obj)
Definition: FLA_Query.c:116
Definition: blis_type_defs.h:137