e99932856ff483f61f274a1abf3d2118e0acc340
[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 #include <math.h>
35 #ifdef __cplusplus
36 extern "C" {
37 #endif
38 #include "cblas.h"
39 #include "f2c.h"
40 #include "blaswrap.h"
41 #include "clapack.h"
42 #ifdef __cplusplus
43 }
44 #endif
46 /*-----------------------------------------------------------------------------
47 * Timing Setup
48 *----------------------------------------------------------------------------*/
49 struct timespec t0,t1;
50 #define tick()  clock_gettime(CLOCK_MONOTONIC, &t0);
51 #define tock() (clock_gettime(CLOCK_MONOTONIC, &t1), \
52                         t1.tv_sec - t0.tv_sec + (t1.tv_nsec - t0.tv_nsec) / 1e9)
54 /* Auxiliary routines prototypes */
55 void print_int_vector( char* desc, integer n, integer* a );
56 void print_real_matrix(char* desc, doublereal *mat, integer n, integer m, integer ldv ); 
57 double compute_inverse_error(doublereal *matrix_A, doublereal *matrix_A_inverse, integer n, integer ld);
58 void print_performance(integer *matrix_sizes, float *secs_runs, int num_test, int num_run);
60 /* Parameters */
61 #define MAX_ALLOWED_ERR (1e-10)
63 /* Main program */
64 int print_data;
65 int main(int argc, char *argv[]) {
66     /* Locals */
67     float secs, secs_lud_inv, *secs_all_runs;
68     int num_test, num_run, i, j;
69     integer n, n_min, n_inc, info, lwork, lda;
70     double max_abs_error;
71     doublereal lwork_opt, *work, *matrix_A, *matrix_LUD;
72     doublereal *vec_dbr_diagr, *vec_dbr_diag, *vec_dbr_diagl;
73     integer    *vec_int_pivot, *vec_int_work;
74     integer    *matrix_sizes;
75     
76     integer const_int6    = 6;
77     doublereal const_dbr1 = 1.;
78     doublereal const_dbr0 = 0.;
79     integer rand_seed[4] = {10,123,278,3579};  /* numbers must lie between 0 and 4095 */
81     fprintf(stdout, "\n\n===== LAPACK example: matrix inversion through LUD =====\n");
82     if(argc == 1) { /* no command line arguments, use default */
83         num_test = 3;
84         n_min    = 1000;
85         n_inc    = 1000;
86         num_run  = 1;
87         print_data = 0;
88     }
89     else if (argc < 6) {
90         printf("Usage: ludinv <num_test> <start_size> <size_inc> <num_run per size> <print data or not (1 or 0)>\n");
91         exit(0);
92     }
93     else {
94         num_test = atoi(argv[1]);
95         n_min    = atoi(argv[2]);
96         n_inc    = atoi(argv[3]);
97         num_run  = atoi(argv[4]);
98         print_data = atoi(argv[5]);
99     }
100     rand_seed[0] = 4675;
101     
102     secs_all_runs = (float *)malloc(sizeof(float)*num_test*num_run);
103     matrix_sizes  = (integer*)malloc(sizeof(integer)*num_test);
104     
105     for(i=0, n=n_min; i<num_test; i++, n+=n_inc) 
106     {
107         lda = n;
108         matrix_sizes[i] = n;
110         printf("\n\nStart matrix inversion by LUD for matrix size %d. \n", (int)n);
111         
112         /* Operation counts of LUD inversion according to LAPACK Working Note 41 */
113         /*operation_counts = 2.*n*n*n - 1.5*n*n + 2.5*n;  */
114         
115         /* Allocate memory for matrix and vectors */
116         matrix_A      = (doublereal *)__malloc_ddr(sizeof(doublereal)*n*n);
117         matrix_LUD    = (doublereal *)__malloc_ddr(sizeof(doublereal)*n*n);
118         vec_int_pivot = (integer    *)__malloc_ddr(sizeof(integer)*n);     
119         vec_int_work  = (integer    *)__malloc_ddr(sizeof(integer)*n);     
120         vec_dbr_diag  = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);  
121         vec_dbr_diagl = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);  
122         vec_dbr_diagr = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);  
123                 
124         j = 0;
125         while(j<num_run) {
126             /* Generate a N by N random matrix */
127             printf("\nStart test run %d out of %d. \n", (int)(j+1), (int)num_run);
128             printf("Generate %d by %d random matrix. \n", (int)n, (int)n);
129             dlatmr_(&n, &n, "S", rand_seed, "N", vec_dbr_diag, &const_int6, &const_dbr1, 
130                     &const_dbr1, "T", "N", vec_dbr_diagl, &const_int6, &const_dbr1, vec_dbr_diagr, 
131                     &const_int6, &const_dbr1, "N", vec_int_pivot, &n, &n, &const_dbr0,
132                     &const_dbr1, "N", matrix_A, &lda, vec_int_work, &info);
133             if(info != 0) {
134                 printf("Random matrix generation has a problem. Returned error value is %d. Start a new run. \n", (int)info);
135                 //goto START_NEW_RUN;
136                 exit(-1);
137             }
138         
139             memcpy((void *)matrix_LUD, (void *)matrix_A, sizeof(doublereal)*n*n);
140             
141             print_real_matrix("Generated random matrix (double real): ", matrix_LUD, n, n, lda);
143             tick();
144             
145             /* Computes LU factorization */
146             dgetrf_(&n, &n, matrix_LUD, &lda, vec_int_pivot, &info);
147             
148             secs_lud_inv = tock();
149             
150             if(info != 0) {
151                 printf("LU factorization has a problem. Returned error value of dgetrf is %d. Start a new run. \n", (int)info);
152                 //goto START_NEW_RUN;
153                 exit(-1);
154             }
155             
156             printf("LU decomposition finished. Time spent on LUD: %f seconds.\n", secs_lud_inv);
157             
158             /* Print details of LU factorization */
159             print_real_matrix( "Details of LU factorization", matrix_LUD, n, n, lda);
160         
161             /* Print pivot indices */
162             print_int_vector( "Pivot indices", n, vec_int_pivot );
163         
164             tick();
165             
166             /* Query and allocate the optimal workspace */
167             lwork = -1;
168             dgetri_(&n, matrix_LUD, &lda, vec_int_pivot, &lwork_opt, &lwork, &info);        
169             lwork = (integer)lwork_opt;
170             work = (doublereal *)__malloc_ddr( lwork*sizeof(doublereal) );
171             
172             /* Compute inversion of matrix based on LU factorization */
173             dgetri_(&n, matrix_LUD, &lda, vec_int_pivot, work, &lwork, &info);    
174             
175             secs = tock();
176             
177             __free_ddr((void *)work);
178             if(info != 0) {
179                 printf("Matrix inversion through LU factorization has a problem.");
180                 printf("Returned error value of dgetri is %d. Start a new run. \n", (int)info);
181                 //goto START_NEW_RUN;
182                 exit(-1);
183             }
184                     
185             secs_lud_inv += secs;
186             secs_all_runs[i*num_run+j] = secs_lud_inv;
187             printf("Matrix inversion finished. Time spent on inversion: %f seconds.\n", secs);
188             printf("Time spent on LUD and inversion: %f seconds.\n", secs_lud_inv);
189            
190             /* Print solution */
191             print_real_matrix( "Inverse of matrix A through LU factorization:", matrix_LUD, n, n, lda);        
192         
193             /* Compute A*inv(A) and precision error: I-A*inv(A) */
194             max_abs_error = compute_inverse_error(matrix_A, matrix_LUD, n, lda);
195             printf("Maximum error of matrix inversion by LUD: %e\n", max_abs_error);
196             
197             if(max_abs_error > MAX_ALLOWED_ERR) {
198                 printf("Error is out of tolerance range.\n");
199                 exit(-1);
200             }
202             j++;
203             
204 START_NEW_RUN:
205             printf("Current run finished. \n");
206         } /* while(j<num_run) */
207         
208         __free_ddr((void *)matrix_A);
209         __free_ddr((void *)matrix_LUD);
210         __free_ddr((void *)vec_int_pivot);
211         __free_ddr((void *)vec_int_work);
212         __free_ddr((void *)vec_dbr_diag);
213         __free_ddr((void *)vec_dbr_diagl);
214         __free_ddr((void *)vec_dbr_diagr);        
215     }  /* for(i=0, n=n_min; i<num_test; i++, n+=n_inc) */
216     
217     print_performance(matrix_sizes, secs_all_runs, num_test, num_run);
218     
219     printf("Passed.\n");
221     exit( 0 );
222
225 /* Auxiliary routine: printing a vector of integers */
226 void print_int_vector( char* desc, integer n, integer* a ) 
228    int j;
230    if (!print_data)
231      return;
232         
233    printf( "\n %s\n", desc );
234    for( j = 0; j < n; j++ ) printf( " %6d", (int)a[j] );
235    printf( "\n" );
238 /* Auxiliary routine: printing a real matrix */
239 void print_real_matrix(char* desc, doublereal *mat, integer m, integer n, integer ldv ) 
241    int i, j;
243    if (!print_data)
244      return;
245      
246    printf( "\n %s\n", desc );
247    for(i=0; i<m; i++) {
248       for(j=0; j<n; j++) {
249          printf( " %10.5f ", mat[i+j*ldv]);
250       }
251       printf( "\n" );
252   }
256 double compute_inverse_error(doublereal *matrix_A, doublereal *matrix_A_inverse, integer n, integer ld)
258    int i, j;
259    double abs_err, relative_err;
260    doublereal *matrix_A_mult_A_inv;
261    double max_err = 0.;
263    matrix_A_mult_A_inv = (doublereal *)malloc(n*n*sizeof(doublereal));
264    
265    cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,1.0,matrix_A,ld,
266                matrix_A_inverse,ld,0.0,matrix_A_mult_A_inv,ld);
267    
268    /* sqrt(sum(err^2))/N, sum(fabs(err))/N, max(fabs(err))  */
269    for(i=0; i<n; i++) {
270       for(j=0; j<n; j++) {
271         if(i==j) {
272           abs_err = fabs(matrix_A_mult_A_inv[i+j*ld]-1.);
273         }
274         else {
275           abs_err = fabs(matrix_A_mult_A_inv[i+j*ld]);
276         }       
277        
278         if(abs_err > max_err) {
279           max_err = abs_err;
280         }
281       }
282    }
283       
284    free(matrix_A_mult_A_inv);
285    return max_err;
289 void print_performance(integer *matrix_sizes, float *secs_runs, int num_test, int num_run)
291     int i, j;
292     float inst, min, max;
293     float tot, avg;
294     
295     printf("Seconds for all matrix sizes:\n");
296     printf("Matrix size   ");
297     printf("    Min           Max        Average  ");
298     printf("\n");
299     
300     for(i=0; i<num_test; i++)
301     {
302         printf("%11d ", (int)matrix_sizes[i]);
303         tot = 0.;
304         max = 0.;
305         min = 1e100;
306         for(j=0; j<num_run; j++) 
307           {
308             inst = secs_runs[i*num_run+j];            
309             tot += inst;
310             
311             if(inst>max) max = inst;
312             if(inst<min) min = inst;
313           }
314         avg = tot/(float)num_run;
315         printf("%12f%12f%12f\n", min, max, avg);
316     }