[Pomp-commits] r384 - pkg/src

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2010-10-18 15:36:48 +0200 (Mon, 18 Oct 2010)
New Revision: 384

Modified:
   pkg/src/pomp_mat.h
   pkg/src/probe_marginal.c
   pkg/src/probe_nlar.c
Log:

- rework pomp_mat.h to make the LAPACK interface a bit more flexible


Modified: pkg/src/pomp_mat.h
===================================================================
--- pkg/src/pomp_mat.h	2010-10-13 19:37:24 UTC (rev 383)
+++ pkg/src/pomp_mat.h	2010-10-18 13:36:48 UTC (rev 384)
@@ -4,12 +4,20 @@
 #include "pomp_internal.h"
 #include <R_ext/Lapack.h>
 
-static R_INLINE void pomp_backsolve (double *a, int m, int n, double *x, int incx) {
+static R_INLINE void pomp_backsolve (double *a, int m, int n, double *x, int incx, 
+				     char *uplo, char *transpose, char *unit) {
   // Level 2 BLAS triangular-matrix solver
   // DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
-  F77_NAME(dtrsv)("U","N","N",&n,a,&m,x,&incx);
+  F77_NAME(dtrsv)(uplo,transpose,unit,&n,a,&m,x,&incx);
 }
 
+static R_INLINE void pomp_qr_x_inv_r (double *a, int m, int n, double *b, int ldb) {
+  // Level 3 BLAS triangular-matrix solver
+  // DTRSM(SIDE,UPLO,TRANS,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+  double alpha = 1.0;
+  F77_NAME(dtrsm)("Right","Upper","No transpose","Non-unit",&m,&n,&alpha,a,&m,b,&ldb);
+}
+
 static R_INLINE void pomp_qr (double *a, int m, int n, int *pivot, double *tau) {
   int info, j, lwork = -1;
   double work1;
@@ -27,20 +35,20 @@
 
 static R_INLINE void pomp_qrqy (double *c, double *a, double *tau, int m, int n, int k, int left, int tp) {
   int lda, info, lwork = -1;
-  char side, trans;
+  char *side, *trans;
   double work1;
  
-  side = (left) ? 'L' : 'R';
+  side = (left) ? "left" : "right";
   lda = (left) ? m : n;
-  trans = (tp) ? 'T' : 'N';
+  trans = (tp) ? "transpose" : "no transpose";
 
   // workspace query
   // DORMQR(SIDE,TRANS,M,N,K,A,LDA,TAU,C,LDC,WORK,LWORK,INFO)
-  F77_NAME(dormqr)(&side,&trans,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info);
+  F77_NAME(dormqr)(side,trans,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info);
   lwork = (int) ceil(work1);
   {				// actual call
     double work[lwork];
-    F77_NAME(dormqr)(&side,&trans,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info);
+    F77_NAME(dormqr)(side,trans,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info);
   }
 }
 

Modified: pkg/src/probe_marginal.c
===================================================================
--- pkg/src/probe_marginal.c	2010-10-13 19:37:24 UTC (rev 383)
+++ pkg/src/probe_marginal.c	2010-10-18 13:36:48 UTC (rev 384)
@@ -136,7 +136,7 @@
 
   // solve R b = Q'x for b
   pomp_qrqy(x,mm,tau,nx,1,np,1,1); // y <- Q'y
-  pomp_backsolve(mm,nx,np,x,1);   // y <- R^{-1} Q'y 
+  pomp_backsolve(mm,nx,np,x,1,"Upper","No transpose","Non-unit");   // y <- R^{-1} y
 
   // unpivot and store the coefficients in beta
   for (i = 0; i < np; i++) beta[pivot[i]] = x[i];

Modified: pkg/src/probe_nlar.c
===================================================================
--- pkg/src/probe_nlar.c	2010-10-13 19:37:24 UTC (rev 383)
+++ pkg/src/probe_nlar.c	2010-10-18 13:36:48 UTC (rev 384)
@@ -121,7 +121,7 @@
       pomp_qr(X,nx,nterms,pivot,tau);
       // then solve R b = Q'y for b
       pomp_qrqy(yp,X,tau,nx,1,nterms,1,1); // y <- Q'y 
-      pomp_backsolve(X,nx,nterms,yp,1);   // y <- R^{-1} Q'y 
+      pomp_backsolve(X,nx,nterms,yp,1,"Upper","No transpose","Non-unit");   // y <- R^{-1} y
       
       // unpivot and store coefficients in beta
       for (i = 0; i < nterms; i++) beta[pivot[i]] = yp[i];