Clean up examples eig and ludinv.
[dense-linear-algebra-libraries/linalg.git] / examples / ludinv / main.c
1 /******************************************************************************
2  * Copyright (c) 2013-2015, Texas Instruments Incorporated - http://www.ti.com/
3  *   All rights reserved.
4  *
5  *   Redistribution and use in source and binary forms, with or without
6  *   modification, are permitted provided that the following conditions are met:
7  *       * Redistributions of source code must retain the above copyright
8  *         notice, this list of conditions and the following disclaimer.
9  *       * Redistributions in binary form must reproduce the above copyright
10  *         notice, this list of conditions and the following disclaimer in the
11  *         documentation and/or other materials provided with the distribution.
12  *       * Neither the name of Texas Instruments Incorporated nor the
13  *         names of its contributors may be used to endorse or promote products
14  *         derived from this software without specific prior written permission.
15  *
16  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17  *   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18  *   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19  *   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
20  *   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21  *   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22  *   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23  *   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24  *   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25  *   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
26  *   THE POSSIBILITY OF SUCH DAMAGE.
27  *****************************************************************************/
29 #include <stdlib.h>
30 #include <stdio.h>
31 #include <string.h>
32 #include <time.h>
33 #ifdef __cplusplus
34 extern "C" {
35 #endif
36 #include "cblas.h"
37 #include "f2c.h"
38 #include "blaswrap.h"
39 #include "clapack.h"
40 #ifdef __cplusplus
41 }
42 #endif
44 /*-----------------------------------------------------------------------------
45 * Timing Setup
46 *----------------------------------------------------------------------------*/
47 struct timespec t0,t1;
48 #define tick()  clock_gettime(CLOCK_MONOTONIC, &t0);
49 #define tock() (clock_gettime(CLOCK_MONOTONIC, &t1), \
50                         t1.tv_sec - t0.tv_sec + (t1.tv_nsec - t0.tv_nsec) / 1e9)
52 /* Auxiliary routines prototypes */
53 void print_int_vector( char* desc, integer n, integer* a );
54 void print_real_matrix(char* desc, doublereal *mat, integer n, integer m, integer ldv ); 
55 double compute_inverse_error(doublereal *matrix_A, doublereal *matrix_A_inverse, integer n, integer ld);
56 void print_performance(integer *matrix_sizes, float *secs_runs, int num_test, int num_run);
58 /* Parameters */
59 #define MAX_ALLOWED_ERR (1e-10)
61 /* Main program */
62 int print_data;
63 int main(int argc, char *argv[]) {
64     /* Locals */
65     float secs, secs_lud_inv, *secs_all_runs;
66     int num_test, num_run, i, j;
67     integer n, n_min, n_inc, info, lwork, lda;
68     double max_abs_error;
69     doublereal lwork_opt, *work, *matrix_A, *matrix_LUD;
70     doublereal *vec_dbr_diagr, *vec_dbr_diag, *vec_dbr_diagl;
71     integer    *vec_int_pivot, *vec_int_work;
72     integer    *matrix_sizes;
73     
74     integer const_int6    = 6;
75     doublereal const_dbr1 = 1.;
76     doublereal const_dbr0 = 0.;
77     integer rand_seed[4] = {10,123,278,3579};  /* numbers must lie between 0 and 4095 */
79     fprintf(stdout, "\n\n===== LAPACK example: matrix inversion through LUD =====\n");
80     if(argc == 1) { /* no command line arguments, use default */
81         num_test = 3;
82         n_min    = 1000;
83         n_inc    = 1000;
84         num_run  = 1;
85         print_data = 0;
86     }
87     else if (argc < 6) {
88         printf("Usage: ludinv <num_test> <start_size> <size_inc> <num_run per size> <print data or not (1 or 0)>\n");
89         exit(0);
90     }
91     else {
92         num_test = atoi(argv[1]);
93         n_min    = atoi(argv[2]);
94         n_inc    = atoi(argv[3]);
95         num_run  = atoi(argv[4]);
96         print_data = atoi(argv[5]);
97     }
98     rand_seed[0] = 4675;
99     
100     secs_all_runs = (float *)malloc(sizeof(float)*num_test*num_run);
101     matrix_sizes  = (integer*)malloc(sizeof(integer)*num_test);
102     
103     for(i=0, n=n_min; i<num_test; i++, n+=n_inc) 
104     {
105         lda = n;
106         matrix_sizes[i] = n;
108         printf("\n\nStart matrix inversion by LUD for matrix size %d. \n", (int)n);
109         
110         /* Operation counts of LUD inversion according to LAPACK Working Note 41 */
111         /*operation_counts = 2.*n*n*n - 1.5*n*n + 2.5*n;  */
112         
113         /* Allocate memory for matrix and vectors */
114         matrix_A      = (doublereal *)__malloc_ddr(sizeof(doublereal)*n*n);
115         matrix_LUD    = (doublereal *)__malloc_ddr(sizeof(doublereal)*n*n);
116         vec_int_pivot = (integer    *)__malloc_ddr(sizeof(integer)*n);     
117         vec_int_work  = (integer    *)__malloc_ddr(sizeof(integer)*n);     
118         vec_dbr_diag  = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);  
119         vec_dbr_diagl = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);  
120         vec_dbr_diagr = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);  
121                 
122         j = 0;
123         while(j<num_run) {
124             /* Generate a N by N random matrix */
125             printf("\nStart test run %d out of %d. \n", (int)(j+1), (int)num_run);
126             printf("Generate %d by %d random matrix. \n", (int)n, (int)n);
127             dlatmr_(&n, &n, "S", rand_seed, "N", vec_dbr_diag, &const_int6, &const_dbr1, 
128                     &const_dbr1, "T", "N", vec_dbr_diagl, &const_int6, &const_dbr1, vec_dbr_diagr, 
129                     &const_int6, &const_dbr1, "N", vec_int_pivot, &n, &n, &const_dbr0,
130                     &const_dbr1, "N", matrix_A, &lda, vec_int_work, &info);
131             if(info != 0) {
132                 printf("Random matrix generation has a problem. Returned error value is %d. Start a new run. \n", (int)info);
133                 //goto START_NEW_RUN;
134                 exit(-1);
135             }
136         
137             memcpy((void *)matrix_LUD, (void *)matrix_A, sizeof(doublereal)*n*n);
138             
139             print_real_matrix("Generated random matrix (double real): ", matrix_LUD, n, n, lda);
141             tick();
142             
143             /* Computes LU factorization */
144             dgetrf_(&n, &n, matrix_LUD, &lda, vec_int_pivot, &info);
145             
146             secs_lud_inv = tock();
147             
148             if(info != 0) {
149                 printf("LU factorization has a problem. Returned error value of dgetrf is %d. Start a new run. \n", (int)info);
150                 //goto START_NEW_RUN;
151                 exit(-1);
152             }
153             
154             printf("LU decomposition finished. Time spent on LUD: %f seconds.\n", (int)secs_lud_inv);
155             
156             /* Print details of LU factorization */
157             print_real_matrix( "Details of LU factorization", matrix_LUD, n, n, lda);
158         
159             /* Print pivot indices */
160             print_int_vector( "Pivot indices", n, vec_int_pivot );
161         
162             tick();
163             
164             /* Query and allocate the optimal workspace */
165             lwork = -1;
166             dgetri_(&n, matrix_LUD, &lda, vec_int_pivot, &lwork_opt, &lwork, &info);        
167             lwork = (integer)lwork_opt;
168             work = (doublereal *)__malloc_ddr( lwork*sizeof(doublereal) );
169             
170             /* Compute inversion of matrix based on LU factorization */
171             dgetri_(&n, matrix_LUD, &lda, vec_int_pivot, work, &lwork, &info);    
172             
173             secs = tock();
174             
175             __free_ddr((void *)work);
176             if(info != 0) {
177                 printf("Matrix inversion through LU factorization has a problem.");
178                 printf("Returned error value of dgetri is %d. Start a new run. \n", (int)info);
179                 //goto START_NEW_RUN;
180                 exit(-1);
181             }
182                     
183             secs_lud_inv += secs;
184             secs_all_runs[i*num_run+j] = secs_lud_inv;
185             printf("Matrix inversion finished. Time spent on inversion: %f seconds.\n", secs);
186             printf("Time spent on LUD and inversion: %f seconds.\n", secs_lud_inv);
187            
188             /* Print solution */
189             print_real_matrix( "Inverse of matrix A through LU factorization:", matrix_LUD, n, n, lda);        
190         
191             /* Compute A*inv(A) and precision error: I-A*inv(A) */
192             max_abs_error = compute_inverse_error(matrix_A, matrix_LUD, n, lda);
193             printf("Maximum error of matrix inversion by LUD: %e\n", max_abs_error);
194             
195             if(max_abs_error > MAX_ALLOWED_ERR) {
196                 printf("Error is out of tolerance range.\n");
197                 exit(-1);
198             }
200             j++;
201             
202 START_NEW_RUN:
203             printf("Current run finished. \n");
204         } /* while(j<num_run) */
205         
206         __free_ddr((void *)matrix_A);
207         __free_ddr((void *)matrix_LUD);
208         __free_ddr((void *)vec_int_pivot);
209         __free_ddr((void *)vec_int_work);
210         __free_ddr((void *)vec_dbr_diag);
211         __free_ddr((void *)vec_dbr_diagl);
212         __free_ddr((void *)vec_dbr_diagr);        
213     }  /* for(i=0, n=n_min; i<num_test; i++, n+=n_inc) */
214     
215     print_performance(matrix_sizes, secs_all_runs, num_test, num_run);
216     
217     printf("Passed.\n");
219     exit( 0 );
220
223 /* Auxiliary routine: printing a vector of integers */
224 void print_int_vector( char* desc, integer n, integer* a ) 
226    int j;
228    if (!print_data)
229      return;
230         
231    printf( "\n %s\n", desc );
232    for( j = 0; j < n; j++ ) printf( " %6d", (int)a[j] );
233    printf( "\n" );
236 /* Auxiliary routine: printing a real matrix */
237 void print_real_matrix(char* desc, doublereal *mat, integer m, integer n, integer ldv ) 
239    int i, j;
241    if (!print_data)
242      return;
243      
244    printf( "\n %s\n", desc );
245    for(i=0; i<m; i++) {
246       for(j=0; j<n; j++) {
247          printf( " %10.5f ", mat[i+j*ldv]);
248       }
249       printf( "\n" );
250   }
254 double compute_inverse_error(doublereal *matrix_A, doublereal *matrix_A_inverse, integer n, integer ld)
256    int i, j;
257    double abs_err, relative_err;
258    doublereal *matrix_A_mult_A_inv;
259    double max_err = 0.;
261    matrix_A_mult_A_inv = (doublereal *)malloc(n*n*sizeof(doublereal));
262    
263    cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,1.0,matrix_A,ld,
264                matrix_A_inverse,ld,0.0,matrix_A_mult_A_inv,ld);
265    
266    /* sqrt(sum(err^2))/N, sum(fabs(err))/N, max(fabs(err))  */
267    for(i=0; i<n; i++) {
268       for(j=0; j<n; j++) {
269         if(i==j) {
270           abs_err = fabs(matrix_A_mult_A_inv[i+j*ld]-1.);
271         }
272         else {
273           abs_err = fabs(matrix_A_mult_A_inv[i+j*ld]);
274         }       
275        
276         if(abs_err > max_err) {
277           max_err = abs_err;
278         }
279       }
280    }
281       
282    free(matrix_A_mult_A_inv);
283    return max_err;
287 void print_performance(integer *matrix_sizes, float *secs_runs, int num_test, int num_run)
289     int i, j;
290     float inst, min, max;
291     float tot, avg;
292     
293     printf("Seconds for all matrix sizes:\n");
294     printf("Matrix size   ");
295     printf("    Min           Max        Average  ");
296     printf("\n");
297     
298     for(i=0; i<num_test; i++)
299     {
300         printf("%11d ", (int)matrix_sizes[i]);
301         tot = 0.;
302         max = 0.;
303         min = 1e100;
304         for(j=0; j<num_run; j++) 
305           {
306             inst = secs_runs[i*num_run+j];            
307             tot += inst;
308             
309             if(inst>max) max = inst;
310             if(inst<min) min = inst;
311           }
312         avg = tot/(float)num_run;
313         printf("%12f%12f%12f\n", min, max, avg);
314     }