Skip to content

Commit

Permalink
introduce FC_LEN_T per recommendation of B. Ripley
Browse files Browse the repository at this point in the history
  • Loading branch information
kingaa committed Sep 26, 2021
1 parent dceb8c5 commit 984501a
Show file tree
Hide file tree
Showing 3 changed files with 16 additions and 11 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
Version: 4.0.4.1
Date: 2021-09-23
Version: 4.0.4.2
Date: 2021-09-24
Authors@R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="[email protected]"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
Expand Down
17 changes: 11 additions & 6 deletions src/pomp_mat.h
Original file line number Diff line number Diff line change
@@ -1,8 +1,13 @@
#ifndef _POMP_MAT_H_
#define _POMP_MAT_H_

#include "pomp_internal.h"
#define USE_FC_LEN_T
#include <Rconfig.h>
#include <R_ext/Lapack.h>
#ifndef FCONE
# define FCONE
#endif
#include "pomp_internal.h"

static R_INLINE void pomp_backsolve (double *a, int lda, int n, double *x, int incx,
char *uplo, char *transpose, char *unit) {
Expand All @@ -14,7 +19,7 @@ static R_INLINE void pomp_backsolve (double *a, int lda, int n, double *x, int i
// UPLO is "U" or "L" depending on whether A is upper or lower triangular
// TRANSPOSE is "T" or "N" depending on whether the transpose is desired
// DIAG is "U" or "N" depending on whether A is unit triangular or not
F77_NAME(dtrsv)(uplo,transpose,unit,&n,a,&lda,x,&incx);
F77_CALL(dtrsv)(uplo,transpose,unit,&n,a,&lda,x,&incx FCONE);
}

static R_INLINE void pomp_qr (double *a, int m, int n, int *pivot, double *tau) {
Expand All @@ -23,11 +28,11 @@ static R_INLINE void pomp_qr (double *a, int m, int n, int *pivot, double *tau)
for (j = 0; j < n; j++) pivot[j] = 0; // zero out the pivots (assumed by DGEQP3)
// LAPACK QR decomposition routine
// DGEQP3(M,N,A,LDA,JPVT,TAU,WORK,LWORK,INFO)
F77_NAME(dgeqp3)(&m,&n,a,&m,pivot,tau,&work1,&lwork,&info); // workspace query
F77_CALL(dgeqp3)(&m,&n,a,&m,pivot,tau,&work1,&lwork,&info); // workspace query
lwork = (int) ceil(work1);
{
double work[lwork];
F77_NAME(dgeqp3)(&m,&n,a,&m,pivot,tau,work,&lwork,&info); // actual call
F77_CALL(dgeqp3)(&m,&n,a,&m,pivot,tau,work,&lwork,&info); // actual call
}
for (j = 0; j < n; j++) pivot[j]--;
}
Expand All @@ -38,11 +43,11 @@ static R_INLINE void pomp_qrqy (double *c, double *a, int lda, double *tau, int

// workspace query
// DORMQR(SIDE,TRANS,M,N,K,A,LDA,TAU,C,LDC,WORK,LWORK,INFO)
F77_NAME(dormqr)(side,transpose,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info);
F77_CALL(dormqr)(side,transpose,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info FCONE);
lwork = (int) ceil(work1);
{ // actual call
double work[lwork];
F77_NAME(dormqr)(side,transpose,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info);
F77_CALL(dormqr)(side,transpose,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info FCONE);
}
}

Expand Down
6 changes: 3 additions & 3 deletions src/synth_lik.c
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,9 @@ static void robust_synth_loglik (double *y, int *dim, double *ydat, double *logl
// do first QR decomposition & backsolve
memcpy(y2,y1,nrow*ncol*sizeof(double));
// LAPACK QR decomposition without pivoting DGEQR2(M,N,A,LDA,TAU,WORK,INFO)
F77_NAME(dgeqr2)(&nrow,&ncol,y2,&nrow,tau,work,&info);
F77_CALL(dgeqr2)(&nrow,&ncol,y2,&nrow,tau,work,&info);
// Level-3 BLAS triangular matrix solver DTRSM(SIDE,UPLO,TRANS,DIAG,M,N,ALPHA,A,LDA,B,LDB)
F77_NAME(dtrsm)("right","upper","no transpose","non-unit",&nrow,&ncol,&one,y2,&nrow,y1,&nrow);
F77_CALL(dtrsm)("right","upper","no transpose","non-unit",&nrow,&ncol,&one,y2,&nrow,y1,&nrow FCONE);

// create Campbell weight vector
d0 = sqrt(ncol)+alpha/sqrt(2.0);
Expand Down Expand Up @@ -90,7 +90,7 @@ static void robust_synth_loglik (double *y, int *dim, double *ydat, double *logl

// do second QR decomposition & backsolve
// LAPACK QR decomposition without pivoting DGEQR2(M,N,A,LDA,TAU,WORK,INFO)
F77_NAME(dgeqr2)(&nrow,&ncol,y1,&nrow,tau,work,&info);
F77_CALL(dgeqr2)(&nrow,&ncol,y1,&nrow,tau,work,&info);
pomp_backsolve(y1,nrow,ncol,ydat,1,"upper","transpose","non-unit");

// compute residual sum of squares and add up logs of diag(R)
Expand Down

0 comments on commit 984501a

Please sign in to comment.