package lapack: enable C interface and add cblas link test

1 files changed, 55 insertions, 0 deletions

diff --git a/src/lapack-test.c b/src/lapack-test.c
new file mode 100644
index 0000000..00bf092
--- /dev/null
+++ b/src/lapack-test.c
@@ -0,0 +1,55 @@
+/*
+ * This file is part of MXE.
+ * See index.html for further information.
+ */
+
+/*  taken from http://www.netlib.org/lapack/lapacke.html */
+
+/*  Calling CGEQRF and CUNGQR to compute Q with workspace querying */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <lapacke_utils.h>
+#include <cblas.h>
+
+int main (int argc, const char * argv[])
+{
+   (void)argc;
+   (void)argv;
+
+   lapack_complex_float *a,*tau,*r,*work,one,zero,query;
+   lapack_int info,m,n,lda,lwork;
+   int i,j;
+   float err;
+   m = 10;   n = 5;   lda = m;
+   one = lapack_make_complex_float(1.0,0.0);
+   zero= lapack_make_complex_float(0.0,0.0);
+   a = calloc(m*n,sizeof(lapack_complex_float));
+   r = calloc(n*n,sizeof(lapack_complex_float));
+   tau = calloc(m,sizeof(lapack_complex_float));
+   for(j=0;j<n;j++)
+      for(i=0;i<m;i++)
+         a[i+j*m] = lapack_make_complex_float(i+1,j+1);
+   info = LAPACKE_cgeqrf_work(LAPACK_COL_MAJOR,m,n,a,lda,tau,&query,-1);
+   lwork = (lapack_int)query;
+   info = LAPACKE_cungqr_work(LAPACK_COL_MAJOR,m,n,n,a,lda,tau,&query,-1);
+   lwork = MAX(lwork,(lapack_int)query);
+   work = calloc(lwork,sizeof(lapack_complex_float));
+   info = LAPACKE_cgeqrf_work(LAPACK_COL_MAJOR,m,n,a,lda,tau,work,lwork);
+   info = LAPACKE_cungqr_work(LAPACK_COL_MAJOR,m,n,n,a,lda,tau,work,lwork);
+   for(j=0;j<n;j++)
+      for(i=0;i<n;i++)
+         r[i+j*n]=(i==j)?-one:zero;
+   cblas_cgemm(CblasColMajor,CblasConjTrans,CblasNoTrans,
+               n,n,m,&one,a,lda,a,lda,&one,r,n);
+   err=0.0;
+   for(i=0;i<n;i++)
+      for(j=0;j<n;j++)
+         err=MAX(err,cabs(r[i+j*n]));
+   printf("error=%e\n",err);
+   free(work);
+   free(tau);
+   free(r);
+   free(a);
+   return(info);
+}