diff options
Diffstat (limited to 'src/lapack-test.c')
-rw-r--r-- | src/lapack-test.c | 55 |
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); +} |