[Bioc-devel] linking LAPACK in windows

[Martin Morgan]

On 04/01/2011 10:07 AM, Cáceres Domínguez, Alejandro wrote:
> Dear list,
>
> I am a new contributor to Bioconductor looking for guidance. My package
> (inveRsion) does not pass the R CMD check in liverpool. Basically, when it
> compiles the vignette, one of the functions returns NaNs when it should
> not. I traced back the bug on a windows machine. Although the package
> installs and compiles the c code without errors, I found that an R
> function ultimately relying on a LAPACK call returns the NaNs (source of
> the bug), when in Mac and Linux it does not. Would this be a problem of
> linking LAPACK in windows?
>
> I've got a Makevars file with the line
>
> PKG_LIBS = $(SUBLIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

Hi Alejandro --

This looks correct to me. When I compile the code below I see

newFreq.c: In function ‘newFreq’:
newFreq.c:36:5: warning: statement with no effect

did you mean

   maxfreq = -1.0;

? Otherwise maxfreq is used uninitialized. I see this also in 
inversionModel.c:94 and when running your test code under valgrind

   R -d valgrind -f test.R

Martin

>
> that seems to work in lamb1 and pelham. I've been back and forth the
> R-extension manual and feel a bit out of my depth. I cannot figure out
> what I should do. I imagine that a way out would be to provide the LAPACK
> routines with the package and compile them with it. My question is: Is
> this an acceptable practice? Otherwise, some light into which sort of
> things I would have to do in the Makevars (Makevars.win?) file will be
> much appreciated. Or am I missing something else?
>
> Please see bellow for particulars.
>
> Many thanks,
> Alejandro
>
> Alejandro Caceres
> Centre de Recerca en Epidemiologia Ambiental
> Doctor Aiguader 88, 0803
> Barcelona, Spain
>
>
> -------------
> newFreq.c compiled with
> -------------
>
> R CMD SHLIB newFreq.c
>
> (with a Makevars in the local dir having the line:
> PKG_LIBS = $(SUBLIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS))
>
> --------
> R code
> --------
>
> Resp<-c(0.93,0.98,1,0.97,1,1,1,1,0.93,1)
> block<-c(4,1,7,2,5,1,5,1,4,1)
> nsub<-10
> lev<-7
> nlev<-7
>
> dyn.load("newFreq.dll")
> # dyn.load("newFreq.so")
>
> .C("newFreq", as.double(Resp), as.double(block), as.double(lev),
> as.integer(nlev),as.integer(nsub), as.double(rep(0,nlev)))[[6]]
>
> ---------
> C code, file: newFreq.c
> ----------
> #include<stdio.h>
> #include<R.h>
> #include<R_ext/Lapack.h>
> #include<R_ext/BLAS.h>
>
>
> void newFreq( double *Resp, double *block, double *lev, int *nlev, int
> *nsub, double *out)
> {
>
>    int level, subject, row, col, inx;
>    double freqLev[*nlev];
>    double bb[*nlev];
>
>    double *AA= (double *)malloc(((*nsub)*(*nsub))*sizeof(*AA));
>
>    double maxfreq;
>
>    /* compute freqLev */
>    for (level=0; level<*nlev; ++level)
>      freqLev[level]=0.0;
>
>
>    for (subject=0; subject<*nsub; ++subject)
>    {
>      for (level=0; level<*nlev; ++level)
>      {
>        if(block[subject]==lev[level])
>        {
>         freqLev[level]=freqLev[level]+Resp[subject];
>        }
>       }
>     }
>
>    /*register max freq*/
>    maxfreq==-1.0;
>    inx=0;
>    for (level=0; level<*nlev; ++level)
>    {
>      if(freqLev[level]>maxfreq)
>      {
>        maxfreq=freqLev[level];
>        inx=level;
>      }
>     }
>
>
>      for (row=0; row<*nlev; ++row)
>      {
>          for (col=0; col<*nlev; ++col)
>          {
>            if (row==col)
>            {
>              AA[row*(*nlev)+col]=1.0;
>            }
>            else
>            {
>              AA[row*(*nlev)+col]=0.0;
>            }
>
>            if (row==inx)
>            {
>              AA[row*(*nlev)+col]=-freqLev[col]/maxfreq;
>            }
>
>            if (col==inx)
>            {
>              AA[row*(*nlev)+col]=1;
>            }
>          }
>
>      }
>
>      /*define bb*/
>      for (row=0; row<*nlev; ++row)
>      {
>        bb[row]=0.0;
>        if(row==inx)
>          bb[inx]=1;
>      }
>
>
>      /*solve linear algebra*/
>      int nr=1;
>      int *pnr;
>      pnr=&nr;
>      *pnr=nr;
>
>      int inf=0;
>      int *pinf;
>      pinf=&inf;
>      *pinf=inf;
>
>      int ipiv[*nlev];
>
>      F77_NAME(dgesv)(nlev, pnr, AA, nlev, ipiv, bb, nlev, pinf);
>
>      for (row=0; row<*nlev; ++row)
>       out[row]=bb[row];
>
>     free(AA);
>     AA=NULL;
>
> }
>
>
>
>
>
>
>
>
>
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> Bioc-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/bioc-devel


-- 
Computational Biology
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109

Location: M1-B861
Telephone: 206 667-2793