rpms
/
R-lpSolve
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
R-lpSolve/0001-Use-R-provided-BLAS-ro...

1023 lines
27 KiB
Diff
Raw Permalink Blame History

From 0acc966a039c1b7646f471b1e75ae37141735ec7 Mon Sep 17 00:00:00 2001
From: Elliott Sales de Andrade <quantum.analyst@gmail.com>
Date: Sun, 27 Sep 2020 05:40:29 -0400
Subject: [PATCH] Use R-provided BLAS routines.
Signed-off-by: Elliott Sales de Andrade <quantum.analyst@gmail.com>
---
src/Makevars | 1 +
src/Makevars.win | 1 +
src/myblas.c | 814 +----------------------------------------------
src/myblas.h | 101 +-----
4 files changed, 12 insertions(+), 905 deletions(-)
diff --git a/src/Makevars b/src/Makevars
index e0c65cc..fb9e2ed 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1 +1,2 @@
PKG_CPPFLAGS=-I . -DINTEGERTIME -DPARSER_LP -DBUILDING_FOR_R -DYY_NEVER_INTERACTIVE -DUSRDLL -DCLOCKTIME -DRoleIsExternalInvEngine -DINVERSE_ACTIVE=INVERSE_LUSOL -DINLINE=static -DParanoia
+PKG_LIBS=$(BLAS_LIBS) $(FLIBS)
diff --git a/src/Makevars.win b/src/Makevars.win
index 4405005..0654c66 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -1 +1,2 @@
PKG_CPPFLAGS=-I . -DINTEGERTIME -DPARSER_LP -DBUILDING_FOR_R -DYY_NEVER_INTERACTIVE -DUSRDLL -DCLOCKTIME -DRoleIsExternalInvEngine -DINVERSE_ACTIVE=INVERSE_LUSOL -DWIN32 -DINLINE=static -DParanoia
+PKG_LIBS=$(BLAS_LIBS) $(FLIBS)
diff --git a/src/myblas.c b/src/myblas.c
index cd6227d..710c3e9 100644
--- a/src/myblas.c
+++ b/src/myblas.c
@@ -1,849 +1,51 @@
-#include <stdlib.h>
-#include <stdio.h>
-/*#include <memory.h>*/
-#include <string.h>
-#include <math.h>
#include "myblas.h"
-#ifdef FORTIFY
-# include "lp_fortify.h"
-#endif
-
-/* ************************************************************************ */
-/* Initialize BLAS interfacing routines */
-/* ************************************************************************ */
-MYBOOL mustinitBLAS = TRUE;
-#if (defined WIN32) || (defined WIN64)
- HINSTANCE hBLAS = NULL;
-#else
- void *hBLAS = NULL;
-#endif
-
-
-/* ************************************************************************ */
-/* Function pointers for external BLAS library (C base 0) */
-/* ************************************************************************ */
-BLAS_dscal_func *BLAS_dscal;
-BLAS_dcopy_func *BLAS_dcopy;
-BLAS_daxpy_func *BLAS_daxpy;
-BLAS_dswap_func *BLAS_dswap;
-BLAS_ddot_func *BLAS_ddot;
-BLAS_idamax_func *BLAS_idamax;
-BLAS_idamin_func *BLAS_idamin;
-BLAS_dload_func *BLAS_dload;
-BLAS_dnormi_func *BLAS_dnormi;
-
-
-/* ************************************************************************ */
-/* Define the BLAS interfacing routines */
-/* ************************************************************************ */
-
void init_BLAS(void)
{
- if(mustinitBLAS) {
- load_BLAS(NULL);
- mustinitBLAS = FALSE;
- }
}
MYBOOL is_nativeBLAS(void)
{
-#ifdef LoadableBlasLib
- return( (MYBOOL) (hBLAS == NULL) );
-#else
- return( TRUE );
-#endif
+ return( FALSE );
}
MYBOOL load_BLAS(char *libname)
{
- MYBOOL result = TRUE;
-
-#ifdef LoadableBlasLib
- if(hBLAS != NULL) {
- my_FreeLibrary(hBLAS);
- }
-#endif
-
- if(libname == NULL) {
- if(!mustinitBLAS && is_nativeBLAS())
- return( FALSE );
- BLAS_dscal = my_dscal;
- BLAS_dcopy = my_dcopy;
- BLAS_daxpy = my_daxpy;
- BLAS_dswap = my_dswap;
- BLAS_ddot = my_ddot;
- BLAS_idamax = my_idamax;
- BLAS_idamin = my_idamin;
- BLAS_dload = my_dload;
- BLAS_dnormi = my_dnormi;
- if(mustinitBLAS)
- mustinitBLAS = FALSE;
- }
- else {
-#ifdef LoadableBlasLib
- #if (defined WIN32) || (defined WIN64)
- char *blasname = libname;
- #else
- /* First standardize UNIX .SO library name format. */
- char blasname[260];
- if(!so_stdname(blasname, libname, 260))
- return( FALSE );
- #endif
- /* Get a handle to the Windows DLL module. */
- hBLAS = my_LoadLibrary(blasname);
-
- /* If the handle is valid, try to get the function addresses. */
- result = (MYBOOL) (hBLAS != NULL);
- if(result) {
- BLAS_dscal = (BLAS_dscal_func *) my_GetProcAddress(hBLAS, BLAS_prec "scal");
- BLAS_dcopy = (BLAS_dcopy_func *) my_GetProcAddress(hBLAS, BLAS_prec "copy");
- BLAS_daxpy = (BLAS_daxpy_func *) my_GetProcAddress(hBLAS, BLAS_prec "axpy");
- BLAS_dswap = (BLAS_dswap_func *) my_GetProcAddress(hBLAS, BLAS_prec "swap");
- BLAS_ddot = (BLAS_ddot_func *) my_GetProcAddress(hBLAS, BLAS_prec "dot");
- BLAS_idamax = (BLAS_idamax_func *) my_GetProcAddress(hBLAS, "i" BLAS_prec "amax");
- BLAS_idamin = (BLAS_idamin_func *) my_GetProcAddress(hBLAS, "i" BLAS_prec "amin");
-#if 0
- BLAS_dload = (BLAS_dload_func *) my_GetProcAddress(hBLAS, BLAS_prec "load");
- BLAS_dnormi = (BLAS_dnormi_func *) my_GetProcAddress(hBLAS, BLAS_prec "normi");
-#endif
- }
-#endif
- /* Do validation */
- if(!result ||
- ((BLAS_dscal == NULL) ||
- (BLAS_dcopy == NULL) ||
- (BLAS_daxpy == NULL) ||
- (BLAS_dswap == NULL) ||
- (BLAS_ddot == NULL) ||
- (BLAS_idamax == NULL) ||
- (BLAS_idamin == NULL) ||
- (BLAS_dload == NULL) ||
- (BLAS_dnormi == NULL))
- ) {
- load_BLAS(NULL);
- result = FALSE;
- }
- }
- return( result );
+ return( TRUE );
}
MYBOOL unload_BLAS(void)
{
- return( load_BLAS(NULL) );
+ return( TRUE );
}
-
-/* ************************************************************************ */
-/* Now define the unoptimized local BLAS functions */
/* ************************************************************************ */
+
void daxpylpsolve( int n, REAL da, REAL *dx, int incx, REAL *dy, int incy)
{
dx++;
dy++;
- BLAS_daxpy( &n, &da, dx, &incx, dy, &incy);
-}
-void BLAS_CALLMODEL my_daxpy( int *_n, REAL *_da, REAL *dx, int *_incx, REAL *dy, int *_incy)
-{
-
-/* constant times a vector plus a vector.
- uses unrolled loops for increments equal to one.
- jack dongarra, linpack, 3/11/78.
- modified 12/3/93, array[1] declarations changed to array[*] */
-
- int i, ix, iy;
-#ifndef DOFASTMATH
- int m, mp1;
-#endif
- register REAL rda;
- REAL da = *_da;
- int n = *_n, incx = *_incx, incy = *_incy;
-
- if (n <= 0) return;
- if (da == 0.0) return;
-
- dx--;
- dy--;
- ix = 1;
- iy = 1;
- if (incx < 0)
- ix = (-n+1)*incx + 1;
- if (incy < 0)
- iy = (-n+1)*incy + 1;
- rda = da;
-
-/* CPU intensive loop; option to do pointer arithmetic */
-#if defined DOFASTMATH
- {
- REAL *xptr, *yptr;
- for (i = 1, xptr = dx + ix, yptr = dy + iy;
- i <= n; i++, xptr += incx, yptr += incy)
- (*yptr) += rda*(*xptr);
- }
-#else
-
- if (incx==1 && incy==1) goto x20;
-
-/* code for unequal increments or equal increments not equal to 1 */
- for (i = 1; i<=n; i++) {
- dy[iy]+= rda*dx[ix];
- ix+= incx;
- iy+= incy;
- }
- return;
-
-/* code for both increments equal to 1 */
-
-/* clean-up loop */
-x20:
- m = n % 4;
- if (m == 0) goto x40;
- for (i = 1; i<=m; i++)
- dy[i]+= rda*dx[i];
- if(n < 4) return;
-x40:
- mp1 = m + 1;
- for (i = mp1; i<=n; i=i+4) {
- dy[i]+= rda*dx[i];
- dy[i + 1]+= rda*dx[i + 1];
- dy[i + 2]+= rda*dx[i + 2];
- dy[i + 3]+= rda*dx[i + 3];
- }
-#endif
+ F77_CALL(daxpy)( &n, &da, dx, &incx, dy, &incy);
}
-
/* ************************************************************************ */
void dcopylpsolve( int n, REAL *dx, int incx, REAL *dy, int incy)
{
dx++;
dy++;
- BLAS_dcopy( &n, dx, &incx, dy, &incy);
-}
-
-void BLAS_CALLMODEL my_dcopy (int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy)
-{
-
-/* copies a vector, x, to a vector, y.
- uses unrolled loops for increments equal to one.
- jack dongarra, linpack, 3/11/78.
- modified 12/3/93, array[1] declarations changed to array[*] */
-
- int i, ix, iy;
-#ifndef DOFASTMATH
- int m, mp1;
-#endif
- int n = *_n, incx = *_incx, incy = *_incy;
-
- if(n<=0)
- return;
-
- dx--;
- dy--;
- ix = 1;
- iy = 1;
- if(incx<0)
- ix = (-n+1)*incx + 1;
- if(incy<0)
- iy = (-n+1)*incy + 1;
-
-
-/* CPU intensive loop; option to do pointer arithmetic */
-#if defined DOFASTMATH
- {
- REAL *xptr, *yptr;
- for (i = 1, xptr = dx + ix, yptr = dy + iy;
- i <= n; i++, xptr += incx, yptr += incy)
- (*yptr) = (*xptr);
- }
-#else
-
- if(incx==1 && incy==1)
- goto x20;
-
-/* code for unequal increments or equal increments not equal to 1 */
-
- for(i = 1; i<=n; i++) {
- dy[iy] = dx[ix];
- ix+= incx;
- iy+= incy;
- }
- return;
-
-/* code for both increments equal to 1 */
-
-/* version with fast machine copy logic (requires memory.h or string.h) */
-x20:
- m = n % 7;
- if (m == 0) goto x40;
- for (i = 1; i<=m; i++)
- dy[i] = dx[i];
- if (n < 7) return;
-
-x40:
- mp1 = m + 1;
- for (i = mp1; i<=n; i=i+7) {
- dy[i] = dx[i];
- dy[i + 1] = dx[i + 1];
- dy[i + 2] = dx[i + 2];
- dy[i + 3] = dx[i + 3];
- dy[i + 4] = dx[i + 4];
- dy[i + 5] = dx[i + 5];
- dy[i + 6] = dx[i + 6];
- }
-#endif
+ F77_CALL(dcopy)( &n, dx, &incx, dy, &incy);
}
-
/* ************************************************************************ */
-
void dscallpsolve (int n, REAL da, REAL *dx, int incx)
{
dx++;
- BLAS_dscal (&n, &da, dx, &incx);
-}
-
-void BLAS_CALLMODEL my_dscal (int *_n, REAL *_da, REAL *dx, int *_incx)
-{
-
-/* Multiply a vector by a constant.
-
- --Input--
- N number of elements in input vector(s)
- DA double precision scale factor
- DX double precision vector with N elements
- INCX storage spacing between elements of DX
-
- --Output--
- DX double precision result (unchanged if N.LE.0)
-
- Replace double precision DX by double precision DA*DX.
- For I = 0 to N-1, replace DX(IX+I*INCX) with DA * DX(IX+I*INCX),
- where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX. */
-
- int i;
-#ifndef DOFASTMATH
- int m, mp1, ix;
-#endif
- register REAL rda;
- REAL da = *_da;
- int n = *_n, incx = *_incx;
-
- if (n <= 0)
- return;
- rda = da;
-
- dx--;
-
-/* Optionally do fast pointer arithmetic */
-#if defined DOFASTMATH
- {
- REAL *xptr;
- for (i = 1, xptr = dx + 1; i <= n; i++, xptr += incx)
- (*xptr) *= rda;
- }
-#else
-
- if (incx == 1)
- goto x20;
-
-/* Code for increment not equal to 1 */
- ix = 1;
- if (incx < 0)
- ix = (-n+1)*incx + 1;
- for(i = 1; i <= n; i++, ix += incx)
- dx[ix] *= rda;
- return;
-
-/* Code for increment equal to 1. */
-/* Clean-up loop so remaining vector length is a multiple of 5. */
-x20:
- m = n % 5;
- if (m == 0) goto x40;
- for( i = 1; i <= m; i++)
- dx[i] *= rda;
- if (n < 5)
- return;
-x40:
- mp1 = m + 1;
- for(i = mp1; i <= n; i += 5) {
- dx[i] *= rda;
- dx[i+1] *= rda;
- dx[i+2] *= rda;
- dx[i+3] *= rda;
- dx[i+4] *= rda;
- }
-#endif
-}
-
-
-/* ************************************************************************ */
-
-REAL ddotlpsolve(int n, REAL *dx, int incx, REAL *dy, int incy)
-{
- dx++;
- dy++;
- return( BLAS_ddot (&n, dx, &incx, dy, &incy) );
-}
-
-REAL BLAS_CALLMODEL my_ddot(int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy)
-{
-
-/* forms the dot product of two vectors.
- uses unrolled loops for increments equal to one.
- jack dongarra, linpack, 3/11/78.
- modified 12/3/93, array[1] declarations changed to array[*] */
-
- register REAL dtemp;
- int i, ix, iy;
-#ifndef DOFASTMATH
- int m, mp1;
-#endif
- int n = *_n, incx = *_incx, incy = *_incy;
-
- dtemp = 0.0;
- if (n<=0)
- return( (REAL) dtemp);
-
- dx--;
- dy--;
- ix = 1;
- iy = 1;
- if (incx<0)
- ix = (-n+1)*incx + 1;
- if (incy<0)
- iy = (-n+1)*incy + 1;
-
-/* CPU intensive loop; option to do pointer arithmetic */
-
-#if defined DOFASTMATH
- {
- REAL *xptr, *yptr;
- for (i = 1, xptr = dx + ix, yptr = dy + iy;
- i <= n; i++, xptr += incx, yptr += incy)
- dtemp+= (*yptr)*(*xptr);
- }
-#else
-
- if (incx==1 && incy==1) goto x20;
-
-/* code for unequal increments or equal increments not equal to 1 */
-
- for (i = 1; i<=n; i++) {
- dtemp+= dx[ix]*dy[iy];
- ix+= incx;
- iy+= incy;
- }
- return(dtemp);
-
-/* code for both increments equal to 1 */
-
-/* clean-up loop */
-
-x20:
- m = n % 5;
- if (m == 0) goto x40;
- for (i = 1; i<=m; i++)
- dtemp+= dx[i]*dy[i];
- if (n < 5) goto x60;
-
-x40:
- mp1 = m + 1;
- for (i = mp1; i<=n; i=i+5)
- dtemp+= dx[i]*dy[i] + dx[i + 1]*dy[i + 1] +
- dx[i + 2]*dy[i + 2] + dx[i + 3]*dy[i + 3] + dx[i + 4]*dy[i + 4];
-
-x60:
-#endif
- return(dtemp);
-}
-
-
-/* ************************************************************************ */
-
-void dswaplpsolve( int n, REAL *dx, int incx, REAL *dy, int incy )
-{
- dx++;
- dy++;
- BLAS_dswap( &n, dx, &incx, dy, &incy );
-}
-
-void BLAS_CALLMODEL my_dswap( int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy )
-{
- int i, ix, iy;
-#ifndef DOFASTMATH
- int m, mp1, ns;
- REAL dtemp2, dtemp3;
-#endif
- register REAL dtemp1;
- int n = *_n, incx = *_incx, incy = *_incy;
-
- if (n <= 0) return;
-
- dx--;
- dy--;
- ix = 1;
- iy = 1;
- if (incx < 0)
- ix = (-n+1)*incx + 1;
- if (incy < 0)
- iy = (-n+1)*incy + 1;
-
-/* CPU intensive loop; option to do pointer arithmetic */
-#if defined DOFASTMATH
- {
- REAL *xptr, *yptr;
- for (i = 1, xptr = dx + ix, yptr = dy + iy;
- i <= n; i++, xptr += incx, yptr += incy) {
- dtemp1 = (*xptr);
- (*xptr) = (*yptr);
- (*yptr) = dtemp1;
- }
- }
-#else
-
- if (incx == incy) {
- if (incx <= 0) goto x5;
- if (incx == 1) goto x20;
- goto x60;
- }
-
-/* code for unequal or nonpositive increments. */
-x5:
- for (i = 1; i<=n; i++) {
- dtemp1 = dx[ix];
- dx[ix] = dy[iy];
- dy[iy] = dtemp1;
- ix+= incx;
- iy+= incy;
- }
- return;
-
-/* code for both increments equal to 1.
- clean-up loop so remaining vector length is a multiple of 3. */
-x20:
- m = n % 3;
- if (m == 0) goto x40;
- for (i = 1; i<=m; i++) {
- dtemp1 = dx[i];
- dx[i] = dy[i];
- dy[i] = dtemp1;
- }
- if (n < 3) return;
-
-x40:
- mp1 = m + 1;
- for (i = mp1; i<=n; i=i+3) {
- dtemp1 = dx[i];
- dtemp2 = dx[i+1];
- dtemp3 = dx[i+2];
- dx[i] = dy[i];
- dx[i+1] = dy[i+1];
- dx[i+2] = dy[i+2];
- dy[i] = dtemp1;
- dy[i+1] = dtemp2;
- dy[i+2] = dtemp3;
- }
- return;
-
-/* code for equal, positive, non-unit increments. */
-x60:
- ns = n*incx;
- for (i = 1; i<=ns; i=i+incx) {
- dtemp1 = dx[i];
- dx[i] = dy[i];
- dy[i] = dtemp1;
- }
-#endif
-
-}
-
-
-/* ************************************************************************ */
-
-void dload(int n, REAL da, REAL *dx, int incx)
-{
- dx++;
- BLAS_dload (&n, &da, dx, &incx);
-}
-
-void BLAS_CALLMODEL my_dload (int *_n, REAL *_da, REAL *dx, int *_incx)
-{
-/* copies a scalar, a, to a vector, x.
- uses unrolled loops when incx equals one.
-
- To change the precision of this program, run the change
- program on dload.f
- Alternatively, to make a single precision version append a
- comment character to the start of all lines between sequential
- precision > double
- and
- end precision > double
- comments and delete the comment character at the start of all
- lines between sequential
- precision > single
- and
- end precision > single
- comments. To make a double precision version interchange
- the append and delete operations in these instructions. */
-
- int i, ix, m, mp1;
- REAL da = *_da;
- int n = *_n, incx = *_incx;
-
- if (n<=0) return;
- dx--;
- if (incx==1) goto x20;
-
-/* code for incx not equal to 1 */
-
- ix = 1;
- if (incx<0)
- ix = (-n+1)*incx + 1;
- for (i = 1; i<=n; i++) {
- dx[ix] = da;
- ix+= incx;
- }
- return;
-
-/* code for incx equal to 1 and clean-up loop */
-
-x20:
- m = n % 7;
- if (m == 0) goto x40;
- for (i = 1; i<=m; i++)
- dx[i] = da;
- if (n < 7) return;
-
-x40:
- mp1 = m + 1;
- for (i = mp1; i<=n; i=i+7) {
- dx[i] = da;
- dx[i + 1] = da;
- dx[i + 2] = da;
- dx[i + 3] = da;
- dx[i + 4] = da;
- dx[i + 5] = da;
- dx[i + 6] = da;
- }
+ F77_CALL(dscal)(&n, &da, dx, &incx);
}
/* ************************************************************************ */
int idamaxlpsolve( int n, REAL *x, int is )
{
x++;
- return ( BLAS_idamax( &n, x, &is ) );
-}
-
-int idaminlpsolve( int n, REAL *x, int is )
-{
- x++;
- return ( BLAS_idamin( &n, x, &is ) );
-}
-
-int BLAS_CALLMODEL my_idamax( int *_n, REAL *x, int *_is )
-{
- register REAL xmax, xtest;
- int i, imax = 0;
-#if !defined DOFASTMATH
- int ii;
-#endif
- int n = *_n, is = *_is;
-
- if((n < 1) || (is <= 0))
- return(imax);
- imax = 1;
- if(n == 1)
- return(imax);
-
-#if defined DOFASTMATH
- xmax = fabs(*x);
- for (i = 2, x += is; i <= n; i++, x += is) {
- xtest = fabs(*x);
- if(xtest > xmax) {
- xmax = xtest;
- imax = i;
- }
- }
-#else
- x--;
- ii = 1;
- xmax = fabs(x[ii]);
- for(i = 2, ii+ = is; i <= n; i++, ii+ = is) {
- xtest = fabs(x[ii]);
- if(xtest > xmax) {
- xmax = xtest;
- imax = i;
- }
- }
-#endif
- return(imax);
-}
-
-int BLAS_CALLMODEL my_idamin( int *_n, REAL *x, int *_is )
-{
- register REAL xmin, xtest;
- int i, imin = 0;
-#if !defined DOFASTMATH
- int ii;
-#endif
- int n = *_n, is = *_is;
-
- if((n < 1) || (is <= 0))
- return(imin);
- imin = 1;
- if(n == 1)
- return(imin);
-
-#if defined DOFASTMATH
- xmin = fabs(*x);
- for (i = 2, x += is; i <= n; i++, x += is) {
- xtest = fabs(*x);
- if(xtest < xmin) {
- xmin = xtest;
- imin = i;
- }
- }
-#else
- x--;
- ii = 1;
- xmin = fabs(x[ii]);
- for(i = 2, ii+ = is; i <= n; i++, ii+ = is) {
- xtest = fabs(x[ii]);
- if(xtest < xmin) {
- xmin = xtest;
- imin = i;
- }
- }
-#endif
- return(imin);
+ return ( F77_CALL(idamax)( &n, x, &is ) );
}
-
-/* ************************************************************************ */
-REAL dnormi( int n, REAL *x )
-{
- x++;
- return( BLAS_dnormi( &n, x ) );
-}
-
-REAL BLAS_CALLMODEL my_dnormi( int *_n, REAL *x )
-{
-/* ===============================================================
- dnormi returns the infinity-norm of the vector x.
- =============================================================== */
- int j;
- register REAL hold, absval;
- int n = *_n;
-
- x--;
- hold = 0.0;
-/* for(j = 1; j <= n; j++) */
- for(j = n; j > 0; j--) {
- absval = fabs(x[j]);
- hold = MAX( hold, absval );
- }
-
- return( hold );
-}
-
-
-/* ************************************************************************ */
-/* Subvector and submatrix access routines (Fortran compatibility) */
-/* ************************************************************************ */
-
-#ifndef UseMacroVector
-int subvec( int item)
-{
- return( item-1 );
-}
-#endif
-
-int submat( int nrowb, int row, int col)
-{
- return( nrowb*(col-1) + subvec(row) );
-}
-
-int posmat( int nrowb, int row, int col)
-{
- return( submat(nrowb, row, col)+BLAS_BASE );
-}
-
-/* ************************************************************************ */
-/* Randomization functions */
-/* ************************************************************************ */
-
-void randomseed(int seeds[])
-/* Simply create some default seed values */
-{
- seeds[1] = 123456;
- seeds[2] = 234567;
- seeds[3] = 345678;
-}
-
-void randomdens( int n, REAL *x, REAL r1, REAL r2, REAL densty, int *seeds )
-{
-/* ------------------------------------------------------------------
- random generates a vector x[*] of random numbers
- in the range (r1, r2) with (approximate) specified density.
- seeds[*] must be initialized before the first call.
- ------------------------------------------------------------------ */
-
- int i;
- REAL *y;
-
- y = (REAL *) malloc(sizeof(*y) * (n+1));
- ddrand( n, x, 1, seeds );
- ddrand( n, y, 1, seeds );
-
- for (i = 1; i<=n; i++) {
- if (y[i] < densty)
- x[i] = r1 + (r2 - r1) * x[i];
- else
- x[i] = 0.0;
- }
- free(y);
-}
-
-
-/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-
-void ddrand( int n, REAL *x, int incx, int *seeds )
-{
-
-/* ------------------------------------------------------------------
- ddrand fills a vector x with uniformly distributed random numbers
- in the interval (0, 1) using a method due to Wichman and Hill.
-
- seeds[1..3] should be set to integer values
- between 1 and 30000 before the first entry.
-
- Integer arithmetic up to 30323 is required.
-
- Blatantly copied from Wichman and Hill 19-January-1987.
- 14-Feb-94. Original version.
- 30 Jun 1999. seeds stored in an array.
- 30 Jun 1999. This version of ddrand.
- ------------------------------------------------------------------ */
-
- int ix, xix;
-
- if (n < 1) return;
-
- for (ix = 1; ix<=1+(n-1)*incx; ix=ix+incx) {
- seeds[1] = 171*(seeds[1] % 177) - 2*(seeds[1]/177);
- seeds[2] = 172*(seeds[2] % 176) - 35*(seeds[2]/176);
- seeds[3] = 170*(seeds[3] % 178) - 63*(seeds[3]/178);
-
- if (seeds[1] < 0) seeds[1] = seeds[1] + 30269;
- if (seeds[2] < 0) seeds[2] = seeds[2] + 30307;
- if (seeds[3] < 0) seeds[3] = seeds[3] + 30323;
-
- x[ix] = ((REAL) seeds[1])/30269.0 +
- ((REAL) seeds[2])/30307.0 +
- ((REAL) seeds[3])/30323.0;
- xix = (int) x[ix];
- x[ix] = fabs(x[ix] - xix);
- }
-
-}
-
diff --git a/src/myblas.h b/src/myblas.h
index 03c68bd..50179e9 100644
--- a/src/myblas.h
+++ b/src/myblas.h
@@ -1,74 +1,18 @@
#ifndef HEADER_myblas
#define HEADER_myblas
-/* ************************************************************************ */
-/* BLAS function interface with local and external loadable versions */
-/* Author: Kjell Eikland */
-/* Version: Initial version spring 2004 */
-/* Licence: LGPL */
-/* ************************************************************************ */
-/* Changes: 19 September 2004 Moved function pointer variable */
-/* declarations from myblas.h to myblas.c */
-/* to avoid linker problems with the Mac. */
-/* 20 April 2005 Modified all double types to REAL to self- */
-/* adjust to global settings. Note that BLAS */
-/* as of now does not have double double. */
-/* 15 December 2005 Added idamin() */
-/* ************************************************************************ */
-
+#include <R_ext/BLAS.h>
+#include "lp_types.h"
#define BLAS_BASE 1
-#define UseMacroVector
-#define LoadableBlasLib
-
-
-/* ************************************************************************ */
-/* Include necessary libraries */
-/* ************************************************************************ */
-#include "commonlib.h"
-#ifdef LoadableBlasLib
- #if (defined WIN32) || (defined WIN64)
- #include <windows.h>
- #else
- #include <dlfcn.h>
- #endif
-#endif
-
#ifdef __cplusplus
extern "C" {
#endif
-
/* ************************************************************************ */
/* BLAS functions */
/* ************************************************************************ */
-#ifndef BLAS_CALLMODEL
- #if (defined WIN32) || (defined WIN64)
- #define BLAS_CALLMODEL _cdecl
- #else
- #define BLAS_CALLMODEL
- #endif
-#endif
-
-typedef void (BLAS_CALLMODEL BLAS_dscal_func) (int *n, REAL *da, REAL *dx, int *incx);
-typedef void (BLAS_CALLMODEL BLAS_dcopy_func) (int *n, REAL *dx, int *incx, REAL *dy, int *incy);
-typedef void (BLAS_CALLMODEL BLAS_daxpy_func) (int *n, REAL *da, REAL *dx, int *incx, REAL *dy, int *incy);
-typedef void (BLAS_CALLMODEL BLAS_dswap_func) (int *n, REAL *dx, int *incx, REAL *dy, int *incy);
-typedef double (BLAS_CALLMODEL BLAS_ddot_func) (int *n, REAL *dx, int *incx, REAL *dy, int *incy);
-typedef int (BLAS_CALLMODEL BLAS_idamax_func)(int *n, REAL *x, int *is);
-typedef int (BLAS_CALLMODEL BLAS_idamin_func)(int *n, REAL *x, int *is);
-typedef void (BLAS_CALLMODEL BLAS_dload_func) (int *n, REAL *da, REAL *dx, int *incx);
-typedef double (BLAS_CALLMODEL BLAS_dnormi_func)(int *n, REAL *x);
-
-#ifndef __WINAPI
- #if (defined WIN32) || (defined WIN64)
- #define __WINAPI WINAPI
- #else
- #define __WINAPI
- #endif
-#endif
-
void init_BLAS(void);
MYBOOL is_nativeBLAS(void);
MYBOOL load_BLAS(char *libname);
@@ -80,48 +24,7 @@ MYBOOL unload_BLAS(void);
void dscallpsolve ( int n, REAL da, REAL *dx, int incx );
void dcopylpsolve ( int n, REAL *dx, int incx, REAL *dy, int incy );
void daxpylpsolve ( int n, REAL da, REAL *dx, int incx, REAL *dy, int incy );
-void dswaplpsolve ( int n, REAL *dx, int incx, REAL *dy, int incy );
-REAL ddotlpsolve ( int n, REAL *dx, int incx, REAL *dy, int incy );
int idamaxlpsolve( int n, REAL *x, int is );
-int idaminlpsolve( int n, REAL *x, int is );
-void dload ( int n, REAL da, REAL *dx, int incx );
-REAL dnormi( int n, REAL *x );
-
-
-/* ************************************************************************ */
-/* Locally implemented BLAS functions (C base 0) */
-/* ************************************************************************ */
-void BLAS_CALLMODEL my_dscal ( int *n, REAL *da, REAL *dx, int *incx );
-void BLAS_CALLMODEL my_dcopy ( int *n, REAL *dx, int *incx, REAL *dy, int *incy );
-void BLAS_CALLMODEL my_daxpy ( int *n, REAL *da, REAL *dx, int *incx, REAL *dy, int *incy );
-void BLAS_CALLMODEL my_dswap ( int *n, REAL *dx, int *incx, REAL *dy, int *incy );
-REAL BLAS_CALLMODEL my_ddot ( int *n, REAL *dx, int *incx, REAL *dy, int *incy );
-int BLAS_CALLMODEL my_idamax( int *n, REAL *x, int *is );
-int BLAS_CALLMODEL my_idamin( int *n, REAL *x, int *is );
-void BLAS_CALLMODEL my_dload ( int *n, REAL *da, REAL *dx, int *incx );
-REAL BLAS_CALLMODEL my_dnormi( int *n, REAL *x );
-
-
-/* ************************************************************************ */
-/* Subvector and submatrix access routines (Fortran compatibility) */
-/* ************************************************************************ */
-#ifdef UseMacroVector
- #define subvec(item) (item - 1)
-#else
- int subvec( int item );
-#endif
-
-int submat( int nrowb, int row, int col );
-int posmat( int nrowb, int row, int col );
-
-
-/* ************************************************************************ */
-/* Randomization functions */
-/* ************************************************************************ */
-void randomseed(int *seeds);
-void randomdens( int n, REAL *x, REAL r1, REAL r2, REAL densty, int *seeds);
-void ddrand( int n, REAL *x, int incx, int *seeds );
-
#ifdef __cplusplus
}
--
2.26.2
Reference in New Issue View Git Blame Copy Permalink