1 /*
2
3 Copyright (C) 2014, The University of Texas at Austin
4
5 This file is part of libflame and is available under the 3-Clause
6 BSD license, which can be found in the LICENSE file at the top-level
7 directory, or at http://opensource.org/licenses/BSD-3-Clause
8
9 */
10
11 #include "FLAME.h"
12 #include "test_libflame.h"
13
14 #define NUM_PARAM_COMBOS 2
15 #define NUM_MATRIX_ARGS 1
16 #define FIRST_VARIANT 1
17 #define LAST_VARIANT 1
18
19 // Static variables.
20 static char* op_str = "Symmetric/Hermitian positive definite inversion";
21 static char* flash_front_str = "FLASH_SPDinv";
22 static char* fla_front_str = "FLA_SPDinv";
23 //static char* fla_unb_var_str = "";
24 //static char* fla_opt_var_str = "";
25 //static char* fla_blk_var_str = "";
26 static char* pc_str[NUM_PARAM_COMBOS] = { "l", "u" };
27 static test_thresh_t thresh = { 1e-02, 1e-03, // warn, pass for s
28 1e-11, 1e-12, // warn, pass for d
29 1e-02, 1e-03, // warn, pass for c
30 1e-11, 1e-12 }; // warn, pass for z
31
32 // Local prototypes.
33 void libfla_test_spdinv_experiment( test_params_t params,
34 unsigned int var,
35 char* sc_str,
36 FLA_Datatype datatype,
37 unsigned int p_cur,
38 unsigned int pci,
39 unsigned int n_repeats,
40 signed int impl,
41 double* perf,
42 double* residual );
43 void libfla_test_spdinv_impl( int impl,
44 FLA_Uplo uplo,
45 FLA_Obj A );
46
47
libfla_test_spdinv(FILE * output_stream,test_params_t params,test_op_t op)48 void libfla_test_spdinv( FILE* output_stream, test_params_t params, test_op_t op )
49 {
50 libfla_test_output_info( "--- %s ---\n", op_str );
51 libfla_test_output_info( "\n" );
52
53 if ( op.flash_front == ENABLE )
54 {
55 //libfla_test_output_info( "%s() front-end...\n", flash_front_str );
56 //libfla_test_output_info( "\n" );
57 libfla_test_op_driver( flash_front_str, NULL,
58 FIRST_VARIANT, LAST_VARIANT,
59 NUM_PARAM_COMBOS, pc_str,
60 NUM_MATRIX_ARGS,
61 FLA_TEST_HIER_FRONT_END,
62 params, thresh, libfla_test_spdinv_experiment );
63 }
64
65 if ( op.fla_front == ENABLE )
66 {
67 //libfla_test_output_info( "%s() front-end...\n", fla_front_str );
68 //libfla_test_output_info( "\n" );
69 libfla_test_op_driver( fla_front_str, NULL,
70 FIRST_VARIANT, LAST_VARIANT,
71 NUM_PARAM_COMBOS, pc_str,
72 NUM_MATRIX_ARGS,
73 FLA_TEST_FLAT_FRONT_END,
74 params, thresh, libfla_test_spdinv_experiment );
75 }
76 }
77
78
79
libfla_test_spdinv_experiment(test_params_t params,unsigned int var,char * sc_str,FLA_Datatype datatype,unsigned int p_cur,unsigned int pci,unsigned int n_repeats,signed int impl,double * perf,double * residual)80 void libfla_test_spdinv_experiment( test_params_t params,
81 unsigned int var,
82 char* sc_str,
83 FLA_Datatype datatype,
84 unsigned int p_cur,
85 unsigned int pci,
86 unsigned int n_repeats,
87 signed int impl,
88 double* perf,
89 double* residual )
90 {
91 dim_t b_flash = params.b_flash;
92 double time_min = 1e9;
93 double time;
94 unsigned int i;
95 unsigned int m;
96 signed int m_input = -1;
97 FLA_Uplo uplo;
98 FLA_Obj A, x, b, norm;
99 FLA_Obj A_save;
100 FLA_Obj A_test, x_test, b_test;
101
102 // Determine the dimensions.
103 if ( m_input < 0 ) m = p_cur / abs(m_input);
104 else m = p_cur;
105
106 // Translate parameter characters to libflame constants.
107 FLA_Param_map_char_to_flame_uplo( &pc_str[pci][0], &uplo );
108
109 // Create the matrices for the current operation.
110 libfla_test_obj_create( datatype, FLA_NO_TRANSPOSE, sc_str[0], m, m, &A );
111
112 // Initialize the test matrices.
113 FLA_Random_spd_matrix( uplo, A );
114
115 // Save the original object contents in a temporary object.
116 FLA_Obj_create_copy_of( FLA_NO_TRANSPOSE, A, &A_save );
117
118 // Create vectors to form a linear system.
119 FLA_Obj_create( datatype, m, 1, 0, 0, &x );
120 FLA_Obj_create( datatype, m, 1, 0, 0, &b );
121
122 // Create a real scalar object to hold the norm of A.
123 FLA_Obj_create( FLA_Obj_datatype_proj_to_real( A ), 1, 1, 0, 0, &norm );
124
125 // Create a random right-hand side vector.
126 FLA_Random_matrix( b );
127
128 // Use hierarchical matrices if we're testing the FLASH front-end.
129 if ( impl == FLA_TEST_HIER_FRONT_END )
130 {
131 FLASH_Obj_create_hier_copy_of_flat( A, 1, &b_flash, &A_test );
132 FLASH_Obj_create_hier_copy_of_flat( b, 1, &b_flash, &b_test );
133 FLASH_Obj_create_hier_copy_of_flat( x, 1, &b_flash, &x_test );
134 }
135 else
136 {
137 A_test = A;
138 b_test = b;
139 x_test = x;
140 }
141
142 // Repeat the experiment n_repeats times and record results.
143 for ( i = 0; i < n_repeats; ++i )
144 {
145 if ( impl == FLA_TEST_HIER_FRONT_END )
146 FLASH_Obj_hierarchify( A_save, A_test );
147 else
148 FLA_Copy_external( A_save, A_test );
149
150 time = FLA_Clock();
151
152 libfla_test_spdinv_impl( impl, uplo, A_test );
153
154 time = FLA_Clock() - time;
155 time_min = min( time_min, time );
156 }
157
158 // Perform a linear solve with the result.
159 if ( impl == FLA_TEST_HIER_FRONT_END )
160 {
161 FLASH_Hemm( FLA_LEFT, uplo,
162 FLA_ONE, A_test, b_test, FLA_ZERO, x_test );
163 FLASH_Obj_flatten( x_test, x );
164 }
165 else
166 {
167 FLA_Hemm( FLA_LEFT, uplo,
168 FLA_ONE, A_test, b_test, FLA_ZERO, x );
169 }
170
171 // Free the hierarchical matrices if we're testing the FLASH front-end.
172 if ( impl == FLA_TEST_HIER_FRONT_END )
173 {
174 FLASH_Obj_free( &A_test );
175 FLASH_Obj_free( &b_test );
176 FLASH_Obj_free( &x_test );
177 }
178
179 // Compute the performance of the best experiment repeat.
180 *perf = 5.0 / 6.0 * m * m * m / time_min / FLOPS_PER_UNIT_PERF;
181 if ( FLA_Obj_is_complex( A ) ) *perf *= 4.0;
182
183 // Compute the residual.
184 FLA_Hemv_external( uplo,
185 FLA_ONE, A_save, x, FLA_MINUS_ONE, b );
186 FLA_Nrm2_external( b, norm );
187 FLA_Obj_extract_real_scalar( norm, residual );
188
189 // Free the supporting flat objects.
190 FLA_Obj_free( &x );
191 FLA_Obj_free( &b );
192 FLA_Obj_free( &norm );
193 FLA_Obj_free( &A_save );
194
195 // Free the flat test matrices.
196 FLA_Obj_free( &A );
197 }
198
199
200
libfla_test_spdinv_impl(int impl,FLA_Uplo uplo,FLA_Obj A)201 void libfla_test_spdinv_impl( int impl,
202 FLA_Uplo uplo,
203 FLA_Obj A )
204 {
205 switch ( impl )
206 {
207 case FLA_TEST_HIER_FRONT_END:
208 FLASH_SPDinv( uplo, A );
209 break;
210
211 case FLA_TEST_FLAT_FRONT_END:
212 FLA_SPDinv( uplo, A );
213 break;
214
215 default:
216 libfla_test_output_error( "Invalid implementation type.\n" );
217 }
218 }
219
220