diff --git a/examples/ludinv/main.c b/examples/ludinv/main.c
index 9dbbd2b99942c1b777e6b98341aaa5a9ef6542ca..106cac15492e444aed6ba85dd13273bac3aadc3a 100644 (file)
--- a/examples/ludinv/main.c
+++ b/examples/ludinv/main.c
extern "C" {
#endif
#include "cblas.h"
+#include "f2c.h"
+#include "blaswrap.h"
+#include "clapack.h"
#ifdef __cplusplus
}
#endif
#define tock() (clock_gettime(CLOCK_MONOTONIC, &t1), \
t1.tv_sec - t0.tv_sec + (t1.tv_nsec - t0.tv_nsec) / 1e9)
-/* LAPACK routines */
-extern int dgetrf_(int *, int *, double *, int *, int *, int *);
-extern int dgetri_(int *, double *, int *, int *, double *, int *, int *);
-extern int dgeqrf_(int *, int *, double *, int *, double *, double *, int *, int *);
-extern int dtrtri_(char *, char *, int *, double *, int *, int *);
-
-extern int dlatmr_(int *, int *, char *, int *, char *,
- double *, int *, double *, double *, char *, char *,
- double *, int *, double *, double *, int *,
- double *, char *, int *, int *, int *,
- double *, double *, char *, double *, int *,
- int *, int *);
-extern int dlacpy_(char *, int *, int *, double *, int *, double *, int *);
-
/* Auxiliary routines prototypes */
-void print_int_vector( char* desc, int n, int* a );
-void print_real_matrix(char* desc, double *mat, int n, int m, int ldv );
-double compute_inverse_error(double *matrix_A, double *matrix_A_inverse, int n, int ld);
-void print_performance(int *matrix_sizes, float *secs_runs, int num_test, int num_run);
+void print_int_vector( char* desc, integer n, integer* a );
+void print_real_matrix(char* desc, doublereal *mat, integer n, integer m, integer ldv );
+double compute_inverse_error(doublereal *matrix_A, doublereal *matrix_A_inverse, integer n, integer ld);
+void print_performance(integer *matrix_sizes, float *secs_runs, int num_test, int num_run);
/* Parameters */
-#define MAX_ALLOWED_ERR (1e-20)
+#define MAX_ALLOWED_ERR (1e-10)
/* Main program */
int print_data;
/* Locals */
float secs, secs_lud_inv, *secs_all_runs;
int num_test, num_run, i, j;
- int n, n_min, n_inc, info, lwork, lda;
- double lwork_opt, max_abs_error;
- double *work;
- double *matrix_A, *matrix_LUD;
- double *vec_dbr_diagr, *vec_dbr_diag, *vec_dbr_diagl;
- int *vec_int_pivot, *vec_int_work;
- int *matrix_sizes;
+ integer n, n_min, n_inc, info, lwork, lda;
+ double max_abs_error;
+ doublereal lwork_opt, *work, *matrix_A, *matrix_LUD;
+ doublereal *vec_dbr_diagr, *vec_dbr_diag, *vec_dbr_diagl;
+ integer *vec_int_pivot, *vec_int_work;
+ integer *matrix_sizes;
- int const_int6 = 6;
- double const_dbr1 = 1.;
- double const_dbr0 = 0.;
- int rand_seed[4] = {10,123,278,3579}; /* numbers must lie between 0 and 4095 */
+ integer const_int6 = 6;
+ doublereal const_dbr1 = 1.;
+ doublereal const_dbr0 = 0.;
+ integer rand_seed[4] = {10,123,278,3579}; /* numbers must lie between 0 and 4095 */
fprintf(stdout, "\n\n===== LAPACK example: matrix inversion through LUD =====\n");
if(argc == 1) { /* no command line arguments, use default */
rand_seed[0] = 4675;
secs_all_runs = (float *)malloc(sizeof(float)*num_test*num_run);
- matrix_sizes = (int*)malloc(sizeof(int)*num_test);
+ matrix_sizes = (integer*)malloc(sizeof(integer)*num_test);
for(i=0, n=n_min; i<num_test; i++, n+=n_inc)
{
lda = n;
matrix_sizes[i] = n;
- printf("\n\nStart matrix inversion by LUD for matrix size %d. \n", n);
+ printf("\n\nStart matrix inversion by LUD for matrix size %d. \n", (int)n);
/* Operation counts of LUD inversion according to LAPACK Working Note 41 */
/*operation_counts = 2.*n*n*n - 1.5*n*n + 2.5*n; */
/* Allocate memory for matrix and vectors */
- matrix_A = (double *)__malloc_ddr(sizeof(double)*n*n);
- matrix_LUD = (double *)__malloc_ddr(sizeof(double)*n*n);
- vec_int_pivot = (int *)__malloc_ddr(sizeof(int)*n);
- vec_int_work = (int *)__malloc_ddr(sizeof(int)*n);
- vec_dbr_diag = (double *)__malloc_ddr(sizeof(double)*n);
- vec_dbr_diagl = (double *)__malloc_ddr(sizeof(double)*n);
- vec_dbr_diagr = (double *)__malloc_ddr(sizeof(double)*n);
+ matrix_A = (doublereal *)__malloc_ddr(sizeof(doublereal)*n*n);
+ matrix_LUD = (doublereal *)__malloc_ddr(sizeof(doublereal)*n*n);
+ vec_int_pivot = (integer *)__malloc_ddr(sizeof(integer)*n);
+ vec_int_work = (integer *)__malloc_ddr(sizeof(integer)*n);
+ vec_dbr_diag = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);
+ vec_dbr_diagl = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);
+ vec_dbr_diagr = (doublereal *)__malloc_ddr(sizeof(doublereal)*n);
j = 0;
while(j<num_run) {
/* Generate a N by N random matrix */
- printf("\nStart test run %d out of %d. \n", j+1, num_run);
- printf("Generate %d by %d random matrix. \n", n, n);
+ printf("\nStart test run %d out of %d. \n", (int)(j+1), (int)num_run);
+ printf("Generate %d by %d random matrix. \n", (int)n, (int)n);
dlatmr_(&n, &n, "S", rand_seed, "N", vec_dbr_diag, &const_int6, &const_dbr1,
&const_dbr1, "T", "N", vec_dbr_diagl, &const_int6, &const_dbr1, vec_dbr_diagr,
&const_int6, &const_dbr1, "N", vec_int_pivot, &n, &n, &const_dbr0,
&const_dbr1, "N", matrix_A, &lda, vec_int_work, &info);
if(info != 0) {
- printf("Random matrix generation has a problem. Returned error value is %d. Start a new run. \n", info);
+ printf("Random matrix generation has a problem. Returned error value is %d. Start a new run. \n", (int)info);
//goto START_NEW_RUN;
exit(-1);
}
- memcpy((void *)matrix_LUD, (void *)matrix_A, sizeof(double)*n*n);
+ memcpy((void *)matrix_LUD, (void *)matrix_A, sizeof(doublereal)*n*n);
print_real_matrix("Generated random matrix (double real): ", matrix_LUD, n, n, lda);
secs_lud_inv = tock();
if(info != 0) {
- printf("LU factorization has a problem. Returned error value of dgetrf is %d. Start a new run. \n", info);
+ printf("LU factorization has a problem. Returned error value of dgetrf is %d. Start a new run. \n", (int)info);
//goto START_NEW_RUN;
exit(-1);
}
- printf("LU decomposition finished. Time spent on LUD: %f seconds.\n", secs_lud_inv);
+ printf("LU decomposition finished. Time spent on LUD: %f seconds.\n", (int)secs_lud_inv);
/* Print details of LU factorization */
print_real_matrix( "Details of LU factorization", matrix_LUD, n, n, lda);
/* Query and allocate the optimal workspace */
lwork = -1;
dgetri_(&n, matrix_LUD, &lda, vec_int_pivot, &lwork_opt, &lwork, &info);
- lwork = (int)lwork_opt;
- work = (double *)__malloc_ddr( lwork*sizeof(double) );
+ lwork = (integer)lwork_opt;
+ work = (doublereal *)__malloc_ddr( lwork*sizeof(doublereal) );
/* Compute inversion of matrix based on LU factorization */
dgetri_(&n, matrix_LUD, &lda, vec_int_pivot, work, &lwork, &info);
__free_ddr((void *)work);
if(info != 0) {
printf("Matrix inversion through LU factorization has a problem.");
- printf("Returned error value of dgetri is %d. Start a new run. \n", info);
+ printf("Returned error value of dgetri is %d. Start a new run. \n", (int)info);
//goto START_NEW_RUN;
exit(-1);
}
/* Compute A*inv(A) and precision error: I-A*inv(A) */
max_abs_error = compute_inverse_error(matrix_A, matrix_LUD, n, lda);
- printf("Maximum relative error of matrix inversion by LUD: %e\n", max_abs_error);
+ printf("Maximum error of matrix inversion by LUD: %e\n", max_abs_error);
if(max_abs_error > MAX_ALLOWED_ERR) {
printf("Error is out of tolerance range.\n");
/* Auxiliary routine: printing a vector of integers */
-void print_int_vector( char* desc, int n, int* a )
+void print_int_vector( char* desc, integer n, integer* a )
{
int j;
return;
printf( "\n %s\n", desc );
- for( j = 0; j < n; j++ ) printf( " %6i", a[j] );
+ for( j = 0; j < n; j++ ) printf( " %6d", (int)a[j] );
printf( "\n" );
}
/* Auxiliary routine: printing a real matrix */
-void print_real_matrix(char* desc, double *mat, int m, int n, int ldv )
+void print_real_matrix(char* desc, doublereal *mat, integer m, integer n, integer ldv )
{
int i, j;
}
-double compute_inverse_error(double *matrix_A, double *matrix_A_inverse, int n, int ld)
+double compute_inverse_error(doublereal *matrix_A, doublereal *matrix_A_inverse, integer n, integer ld)
{
int i, j;
double abs_err, relative_err;
- double *matrix_A_mult_A_inv;
+ doublereal *matrix_A_mult_A_inv;
double max_err = 0.;
- matrix_A_mult_A_inv = (double *)malloc(n*n*sizeof(double));
+ matrix_A_mult_A_inv = (doublereal *)malloc(n*n*sizeof(doublereal));
cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,1.0,matrix_A,ld,
matrix_A_inverse,ld,0.0,matrix_A_mult_A_inv,ld);
- /* sqrt(sum(err^2))/N, sum(abs(err))/N, max(abs(err)) */
+ /* sqrt(sum(err^2))/N, sum(fabs(err))/N, max(fabs(err)) */
for(i=0; i<n; i++) {
for(j=0; j<n; j++) {
if(i==j) {
- abs_err = abs(matrix_A_mult_A_inv[i+j*ld]-1.);
+ abs_err = fabs(matrix_A_mult_A_inv[i+j*ld]-1.);
}
else {
- abs_err = abs(matrix_A_mult_A_inv[i+j*ld]);
- }
-
- relative_err = abs_err / abs(matrix_A[i+j*ld]);
-
- if(relative_err > max_err) {
- max_err = relative_err;
+ abs_err = fabs(matrix_A_mult_A_inv[i+j*ld]);
+ }
+
+ if(abs_err > max_err) {
+ max_err = abs_err;
}
}
}
}
-void print_performance(int *matrix_sizes, float *secs_runs, int num_test, int num_run)
+void print_performance(integer *matrix_sizes, float *secs_runs, int num_test, int num_run)
{
int i, j;
float inst, min, max;
@@ -311,19 +297,19 @@ void print_performance(int *matrix_sizes, float *secs_runs, int num_test, int nu
for(i=0; i<num_test; i++)
{
- printf("%11d ", matrix_sizes[i]);
+ printf("%11d ", (int)matrix_sizes[i]);
tot = 0.;
max = 0.;
min = 1e100;
for(j=0; j<num_run; j++)
{
inst = secs_runs[i*num_run+j];
- tot += (double)inst;
+ tot += inst;
if(inst>max) max = inst;
if(inst<min) min = inst;
}
- avg = tot/(double)num_run;
+ avg = tot/(float)num_run;
printf("%12f%12f%12f\n", min, max, avg);
}
}