diff --git a/com.oracle.truffle.r.native/Makefile b/com.oracle.truffle.r.native/Makefile
index fcc82719fd07527dbca2faa3463ce892557b817d..5ef07e48d29f61366c4d6baa12640d644b92e8be 100644
--- a/com.oracle.truffle.r.native/Makefile
+++ b/com.oracle.truffle.r.native/Makefile
@@ -21,16 +21,38 @@
# questions.
#
-all:
ifneq ($(shell uname), Darwin)
- gcc -fPIC -shared -o ./lib/linux/libRDerived.so ./src/fft.c
+OS_DIR=linux
else
- gcc -fPIC -dynamiclib -o ./lib/darwin/libRDerived.dylib ./src/fft.c
+OS_DIR=darwin
endif
-clean:
+SRC = src
+OBJ = lib/$(OS_DIR)
+
+C_SOURCES := $(wildcard $(SRC)/*.c)
+F_SOURCES := $(wildcard $(SRC)/*.f)
+
+C_OBJECTS := $(subst $(SRC),$(OBJ),$(C_SOURCES:.c=.o))
+F_OBJECTS := $(subst $(SRC),$(OBJ),$(F_SOURCES:.f=.o))
+
+all: $(C_OBJECTS) $(F_OBJECTS)
ifneq ($(shell uname), Darwin)
- rm -f ./lib/linux/libRDerived.*
+ gcc -fPIC -shared -o ./lib/linux/libRDerived.so $(C_OBJECTS) $(F_OBJECTS)
else
- rm -f ./lib/darwin/libRDerived.*
+ gcc -dynamiclib -undefined dynamic_lookup -o ./lib/darwin/libRDerived.dylib $(C_OBJECTS) $(F_OBJECTS)
endif
+
+$(OBJ)/%.o: $(SRC)/%.c
+ gcc -fPIC -O2 -c $< -o $@
+
+$(OBJ)/%.o: $(SRC)/%.f
+ gfortran -fPIC -O2 -c $< -o $@
+
+
+clean:
+ rm -f ./lib/$(OS_DIR)/libRDerived.*
+ rm -f ./lib/$(OS_DIR)/*.o
+
+cleanobj:
+ rm -f ./lib/$(OS_DIR)/*.o
diff --git a/com.oracle.truffle.r.native/lib/darwin/libR.dylib b/com.oracle.truffle.r.native/lib/darwin/libR.dylib
deleted file mode 100755
index 6dea2826530ec11fb3ebd03d113c7fa306b85550..0000000000000000000000000000000000000000
Binary files a/com.oracle.truffle.r.native/lib/darwin/libR.dylib and /dev/null differ
diff --git a/com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib b/com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib
index 8d58a619d4c489fe40556023a711d7e24dc88d08..66e95724987b1332831dd1ae19e6c627c00e511e 100755
Binary files a/com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib and b/com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib differ
diff --git a/com.oracle.truffle.r.native/lib/linux/libRDerived.so b/com.oracle.truffle.r.native/lib/linux/libRDerived.so
index df669c1619b86fcc6e372d68e839821f88c0205d..f1baa25c44c6dfd1638c2092596acd382bf35ba7 100755
Binary files a/com.oracle.truffle.r.native/lib/linux/libRDerived.so and b/com.oracle.truffle.r.native/lib/linux/libRDerived.so differ
diff --git a/com.oracle.truffle.r.native/src/dchdc.f b/com.oracle.truffle.r.native/src/dchdc.f
new file mode 100644
index 0000000000000000000000000000000000000000..81b6786bde7f53ae804c641e6f7a66ce12734ef5
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dchdc.f
@@ -0,0 +1,234 @@
+ subroutine dchdc(a,lda,p,work,jpvt,job,info)
+ integer lda,p,jpvt(p),job,info
+ double precision a(lda,p),work(*)
+c
+c dchdc computes the cholesky decomposition of a positive definite
+c matrix. a pivoting option allows the user to estimate the
+c condition of a positive definite matrix or determine the rank
+c of a positive semidefinite matrix.
+c
+c on entry
+c
+c a double precision(lda,p).
+c a contains the matrix whose decomposition is to
+c be computed. onlt the upper half of a need be stored.
+c the lower part of the array a is not referenced.
+c
+c lda integer.
+c lda is the leading dimension of the array a.
+c
+c p integer.
+c p is the order of the matrix.
+c
+c work double precision.
+c work is a work array.
+c
+c jpvt integer(p).
+c jpvt contains integers that control the selection
+c of the pivot elements, if pivoting has been requested.
+c each diagonal element a(k,k)
+c is placed in one of three classes according to the
+c value of jpvt(k).
+c
+c if jpvt(k) .gt. 0, then x(k) is an initial
+c element.
+c
+c if jpvt(k) .eq. 0, then x(k) is a free element.
+c
+c if jpvt(k) .lt. 0, then x(k) is a final element.
+c
+c before the decomposition is computed, initial elements
+c are moved by symmetric row and column interchanges to
+c the beginning of the array a and final
+c elements to the end. both initial and final elements
+c are frozen in place during the computation and only
+c free elements are moved. at the k-th stage of the
+c reduction, if a(k,k) is occupied by a free element
+c it is interchanged with the largest free element
+c a(l,l) with l .ge. k. jpvt is not referenced if
+c job .eq. 0.
+c
+c job integer.
+c job is an integer that initiates column pivoting.
+c if job .eq. 0, no pivoting is done.
+c if job .ne. 0, pivoting is done.
+c
+c on return
+c
+c a a contains in its upper half the cholesky factor
+c of the matrix a as it has been permuted by pivoting.
+c
+c jpvt jpvt(j) contains the index of the diagonal element
+c of a that was moved into the j-th position,
+c provided pivoting was requested.
+c
+c info contains the index of the last positive diagonal
+c element of the cholesky factor.
+c
+c for positive definite matrices info = p is the normal return.
+c for pivoting with positive semidefinite matrices info will
+c in general be less than p. however, info may be greater than
+c the rank of a, since rounding error can cause an otherwise zero
+c element to be positive. indefinite systems will always cause
+c info to be less than p.
+c
+c linpack. this version dated 08/14/78 .
+c j.j. dongarra and g.w. stewart, argonne national laboratory and
+c university of maryland.
+c
+c
+c blas daxpy,dswap
+c fortran dsqrt
+c
+c internal variables
+c
+ integer pu,pl,plp1,j,jp,jt,k,kb,km1,kp1,l,maxl
+ double precision temp
+ double precision maxdia
+ logical swapk,negk
+c
+ pl = 1
+ pu = 0
+ info = p
+ if (job .eq. 0) go to 160
+c
+c pivoting has been requested. rearrange the
+c the elements according to jpvt.
+c
+ do 70 k = 1, p
+ swapk = jpvt(k) .gt. 0
+ negk = jpvt(k) .lt. 0
+ jpvt(k) = k
+ if (negk) jpvt(k) = -jpvt(k)
+ if (.not.swapk) go to 60
+ if (k .eq. pl) go to 50
+ call dswap(pl-1,a(1,k),1,a(1,pl),1)
+ temp = a(k,k)
+ a(k,k) = a(pl,pl)
+ a(pl,pl) = temp
+ plp1 = pl + 1
+ if (p .lt. plp1) go to 40
+ do 30 j = plp1, p
+ if (j .ge. k) go to 10
+ temp = a(pl,j)
+ a(pl,j) = a(j,k)
+ a(j,k) = temp
+ go to 20
+ 10 continue
+ if (j .eq. k) go to 20
+ temp = a(k,j)
+ a(k,j) = a(pl,j)
+ a(pl,j) = temp
+ 20 continue
+ 30 continue
+ 40 continue
+ jpvt(k) = jpvt(pl)
+ jpvt(pl) = k
+ 50 continue
+ pl = pl + 1
+ 60 continue
+ 70 continue
+ pu = p
+ if (p .lt. pl) go to 150
+ do 140 kb = pl, p
+ k = p - kb + pl
+ if (jpvt(k) .ge. 0) go to 130
+ jpvt(k) = -jpvt(k)
+ if (pu .eq. k) go to 120
+ call dswap(k-1,a(1,k),1,a(1,pu),1)
+ temp = a(k,k)
+ a(k,k) = a(pu,pu)
+ a(pu,pu) = temp
+ kp1 = k + 1
+ if (p .lt. kp1) go to 110
+ do 100 j = kp1, p
+ if (j .ge. pu) go to 80
+ temp = a(k,j)
+ a(k,j) = a(j,pu)
+ a(j,pu) = temp
+ go to 90
+ 80 continue
+ if (j .eq. pu) go to 90
+ temp = a(k,j)
+ a(k,j) = a(pu,j)
+ a(pu,j) = temp
+ 90 continue
+ 100 continue
+ 110 continue
+ jt = jpvt(k)
+ jpvt(k) = jpvt(pu)
+ jpvt(pu) = jt
+ 120 continue
+ pu = pu - 1
+ 130 continue
+ 140 continue
+ 150 continue
+ 160 continue
+ do 270 k = 1, p
+c
+c reduction loop.
+c
+ maxdia = a(k,k)
+ kp1 = k + 1
+ maxl = k
+c
+c determine the pivot element.
+c
+ if (k .lt. pl .or. k .ge. pu) go to 190
+ do 180 l = kp1, pu
+ if (a(l,l) .le. maxdia) go to 170
+ maxdia = a(l,l)
+ maxl = l
+ 170 continue
+ 180 continue
+ 190 continue
+c
+c quit if the pivot element is not positive.
+c
+ if (maxdia .gt. 0.0d0) go to 200
+ info = k - 1
+c ......exit
+ go to 280
+ 200 continue
+ if (k .eq. maxl) go to 210
+c
+c start the pivoting and update jpvt.
+c
+ km1 = k - 1
+ call dswap(km1,a(1,k),1,a(1,maxl),1)
+ a(maxl,maxl) = a(k,k)
+ a(k,k) = maxdia
+ jp = jpvt(maxl)
+ jpvt(maxl) = jpvt(k)
+ jpvt(k) = jp
+ 210 continue
+c
+c reduction step. pivoting is contained across the rows.
+c
+ work(k) = dsqrt(a(k,k))
+ a(k,k) = work(k)
+ if (p .lt. kp1) go to 260
+ do 250 j = kp1, p
+ if (k .eq. maxl) go to 240
+ if (j .ge. maxl) go to 220
+ temp = a(k,j)
+ a(k,j) = a(j,maxl)
+ a(j,maxl) = temp
+ go to 230
+ 220 continue
+ if (j .eq. maxl) go to 230
+ temp = a(k,j)
+ a(k,j) = a(maxl,j)
+ a(maxl,j) = temp
+ 230 continue
+ 240 continue
+ a(k,j) = a(k,j)/work(k)
+ work(j) = a(k,j)
+ temp = -a(k,j)
+ call daxpy(j-k,temp,work(kp1),1,a(kp1,j),1)
+ 250 continue
+ 260 continue
+ 270 continue
+ 280 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dpbfa.f b/com.oracle.truffle.r.native/src/dpbfa.f
new file mode 100644
index 0000000000000000000000000000000000000000..3c22a86094bd3ddb31493f63f09115d0106124fd
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dpbfa.f
@@ -0,0 +1,100 @@
+ subroutine dpbfa(abd,lda,n,m,info)
+
+ integer lda,n,m,info
+ double precision abd(lda,n)
+c
+c dpbfa factors a double precision symmetric positive definite
+c matrix stored in band form.
+c
+c dpbfa is usually called by dpbco, but it can be called
+c directly with a saving in time if rcond is not needed.
+c
+c on entry
+c
+c abd double precision(lda, n)
+c the matrix to be factored. the columns of the upper
+c triangle are stored in the columns of abd and the
+c diagonals of the upper triangle are stored in the
+c rows of abd . see the comments below for details.
+c
+c lda integer
+c the leading dimension of the array abd .
+c lda must be .ge. m + 1 .
+c
+c n integer
+c the order of the matrix a .
+c
+c m integer
+c the number of diagonals above the main diagonal.
+c 0 .le. m .lt. n .
+c
+c on return
+c
+c abd an upper triangular matrix r , stored in band
+c form, so that a = trans(r)*r .
+c
+c info integer
+c = 0 for normal return.
+c = k if the leading minor of order k is not
+c positive definite.
+c
+c band storage
+c
+c if a is a symmetric positive definite band matrix,
+c the following program segment will set up the input.
+c
+c m = (band width above diagonal)
+c do 20 j = 1, n
+c i1 = max0(1, j-m)
+c do 10 i = i1, j
+c k = i-j+m+1
+c abd(k,j) = a(i,j)
+c 10 continue
+c 20 continue
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas ddot
+c fortran max0,sqrt
+c
+c internal variables
+c
+ double precision ddot,t
+ double precision s
+ integer ik,j,jk,k,mu
+c begin block with ...exits to 40
+c
+c
+ do 30 j = 1, n
+ info = j
+ s = 0.0d0
+ ik = m + 1
+ jk = max0(j-m,1)
+ mu = max0(m+2-j,1)
+ if (m .lt. mu) go to 20
+ do 10 k = mu, m
+
+ t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1)
+ t = t/abd(m+1,jk)
+ abd(k,j) = t
+ s = s + t*t
+ ik = ik - 1
+ jk = jk + 1
+ 10 continue
+ 20 continue
+
+ s = abd(m+1,j) - s
+c ......exit
+ if (s .le. 0.0d0) go to 40
+
+ abd(m+1,j) = sqrt(s)
+
+ 30 continue
+ info = 0
+ 40 continue
+ return
+ end
+
diff --git a/com.oracle.truffle.r.native/src/dpbsl.f b/com.oracle.truffle.r.native/src/dpbsl.f
new file mode 100644
index 0000000000000000000000000000000000000000..d910deef88ecd031b13d86af743201e31b49be0d
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dpbsl.f
@@ -0,0 +1,83 @@
+ subroutine dpbsl(abd,lda,n,m,b)
+
+ integer lda,n,m
+ double precision abd(lda,n),b(n)
+c
+c dpbsl solves the double precision symmetric positive definite
+c band system a*x = b
+c using the factors computed by dpbco or dpbfa.
+c
+c on entry
+c
+c abd double precision(lda, n)
+c the output from dpbco or dpbfa.
+c
+c lda integer
+c the leading dimension of the array abd .
+c
+c n integer
+c the order of the matrix a .
+c
+c m integer
+c the number of diagonals above the main diagonal.
+c
+c b double precision(n)
+c the right hand side vector.
+c
+c on return
+c
+c b the solution vector x .
+c
+c error condition
+c
+c a division by zero will occur if the input factor contains
+c a zero on the diagonal. technically this indicates
+c singularity but it is usually caused by improper subroutine
+c arguments. it will not occur if the subroutines are called
+c correctly and info .eq. 0 .
+c
+c to compute inverse(a) * c where c is a matrix
+c with p columns
+c call dpbco(abd,lda,n,rcond,z,info)
+c if (rcond is too small .or. info .ne. 0) go to ...
+c do 10 j = 1, p
+c call dpbsl(abd,lda,n,c(1,j))
+c 10 continue
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas daxpy,ddot
+c fortran min0
+c
+c internal variables
+c
+ double precision ddot,t
+ integer k,kb,la,lb,lm
+c
+c solve trans(r)*y = b
+c
+ do 10 k = 1, n
+ lm = min0(k-1,m)
+ la = m + 1 - lm
+ lb = k - lm
+ t = ddot(lm,abd(la,k),1,b(lb),1)
+ b(k) = (b(k) - t)/abd(m+1,k)
+ 10 continue
+c
+c solve r*x = y
+c
+ do 20 kb = 1, n
+ k = n + 1 - kb
+ lm = min0(k-1,m)
+ la = m + 1 - lm
+ lb = k - lm
+ b(k) = b(k)/abd(m+1,k)
+ t = -b(k)
+ call daxpy(lm,t,abd(la,k),1,b(lb),1)
+ 20 continue
+ return
+ end
+
diff --git a/com.oracle.truffle.r.native/src/dpoco.f b/com.oracle.truffle.r.native/src/dpoco.f
new file mode 100644
index 0000000000000000000000000000000000000000..83a7f1bd9b5f259f6e04549adec7d692a8011200
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dpoco.f
@@ -0,0 +1,195 @@
+c
+c dpoco factors a double precision symmetric positive definite
+c matrix and estimates the condition of the matrix.
+c
+c if rcond is not needed, dpofa is slightly faster.
+c to solve a*x = b , follow dpoco by dposl.
+c to compute inverse(a)*c , follow dpoco by dposl.
+c to compute determinant(a) , follow dpoco by dpodi.
+c to compute inverse(a) , follow dpoco by dpodi.
+c
+c on entry
+c
+c a double precision(lda, n)
+c the symmetric matrix to be factored. only the
+c diagonal and upper triangle are used.
+c
+c lda integer
+c the leading dimension of the array a .
+c
+c n integer
+c the order of the matrix a .
+c
+c on return
+c
+c a an upper triangular matrix r so that a = trans(r)*r
+c where trans(r) is the transpose.
+c the strict lower triangle is unaltered.
+c if info .ne. 0 , the factorization is not complete.
+c
+c rcond double precision
+c an estimate of the reciprocal condition of a .
+c for the system a*x = b , relative perturbations
+c in a and b of size epsilon may cause
+c relative perturbations in x of size epsilon/rcond .
+c if rcond is so small that the logical expression
+c 1.0 + rcond .eq. 1.0
+c is true, then a may be singular to working
+c precision. in particular, rcond is zero if
+c exact singularity is detected or the estimate
+c underflows. if info .ne. 0 , rcond is unchanged.
+c
+c z double precision(n)
+c a work vector whose contents are usually unimportant.
+c if a is close to a singular matrix, then z is
+c an approximate null vector in the sense that
+c norm(a*z) = rcond*norm(a)*norm(z) .
+c if info .ne. 0 , z is unchanged.
+c
+c info integer
+c = 0 for normal return.
+c = k signals an error condition. the leading minor
+c of order k is not positive definite.
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c linpack dpofa
+c blas daxpy,ddot,dscal,dasum
+c fortran dabs,dmax1,dreal,dsign
+c
+ subroutine dpoco(a,lda,n,rcond,z,info)
+ integer lda,n,info
+ double precision a(lda,*),z(*)
+ double precision rcond
+c
+c internal variables
+c
+ double precision ddot,ek,t,wk,wkm
+ double precision anorm,s,dasum,sm,ynorm
+ integer i,j,jm1,k,kb,kp1
+c
+c
+c find norm of a using only upper half
+c
+ do 30 j = 1, n
+ z(j) = dasum(j,a(1,j),1)
+ jm1 = j - 1
+ if (jm1 .lt. 1) go to 20
+ do 10 i = 1, jm1
+ z(i) = z(i) + dabs(a(i,j))
+ 10 continue
+ 20 continue
+ 30 continue
+ anorm = 0.0d0
+ do 40 j = 1, n
+ anorm = dmax1(anorm,z(j))
+ 40 continue
+c
+c factor
+c
+ call dpofa(a,lda,n,info)
+ if (info .ne. 0) go to 180
+c
+c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
+c estimate = norm(z)/norm(y) where a*z = y and a*y = e .
+c the components of e are chosen to cause maximum local
+c growth in the elements of w where trans(r)*w = e .
+c the vectors are frequently rescaled to avoid overflow.
+c
+c solve trans(r)*w = e
+c
+ ek = 1.0d0
+ do 50 j = 1, n
+ z(j) = 0.0d0
+ 50 continue
+ do 110 k = 1, n
+ if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
+ if (dabs(ek-z(k)) .le. a(k,k)) go to 60
+ s = a(k,k)/dabs(ek-z(k))
+ call dscal(n,s,z,1)
+ ek = s*ek
+ 60 continue
+ wk = ek - z(k)
+ wkm = -ek - z(k)
+ s = dabs(wk)
+ sm = dabs(wkm)
+ wk = wk/a(k,k)
+ wkm = wkm/a(k,k)
+ kp1 = k + 1
+ if (kp1 .gt. n) go to 100
+ do 70 j = kp1, n
+ sm = sm + dabs(z(j)+wkm*a(k,j))
+ z(j) = z(j) + wk*a(k,j)
+ s = s + dabs(z(j))
+ 70 continue
+ if (s .ge. sm) go to 90
+ t = wkm - wk
+ wk = wkm
+ do 80 j = kp1, n
+ z(j) = z(j) + t*a(k,j)
+ 80 continue
+ 90 continue
+ 100 continue
+ z(k) = wk
+ 110 continue
+ s = 1.0d0/dasum(n,z,1)
+ call dscal(n,s,z,1)
+c
+c solve r*y = w
+c
+ do 130 kb = 1, n
+ k = n + 1 - kb
+ if (dabs(z(k)) .le. a(k,k)) go to 120
+ s = a(k,k)/dabs(z(k))
+ call dscal(n,s,z,1)
+ 120 continue
+ z(k) = z(k)/a(k,k)
+ t = -z(k)
+ call daxpy(k-1,t,a(1,k),1,z(1),1)
+ 130 continue
+ s = 1.0d0/dasum(n,z,1)
+ call dscal(n,s,z,1)
+c
+ ynorm = 1.0d0
+c
+c solve trans(r)*v = y
+c
+ do 150 k = 1, n
+ z(k) = z(k) - ddot(k-1,a(1,k),1,z(1),1)
+ if (dabs(z(k)) .le. a(k,k)) go to 140
+ s = a(k,k)/dabs(z(k))
+ call dscal(n,s,z,1)
+ ynorm = s*ynorm
+ 140 continue
+ z(k) = z(k)/a(k,k)
+ 150 continue
+ s = 1.0d0/dasum(n,z,1)
+ call dscal(n,s,z,1)
+ ynorm = s*ynorm
+c
+c solve r*z = v
+c
+ do 170 kb = 1, n
+ k = n + 1 - kb
+ if (dabs(z(k)) .le. a(k,k)) go to 160
+ s = a(k,k)/dabs(z(k))
+ call dscal(n,s,z,1)
+ ynorm = s*ynorm
+ 160 continue
+ z(k) = z(k)/a(k,k)
+ t = -z(k)
+ call daxpy(k-1,t,a(1,k),1,z(1),1)
+ 170 continue
+c make znorm = 1.0
+ s = 1.0d0/dasum(n,z,1)
+ call dscal(n,s,z,1)
+ ynorm = s*ynorm
+c
+ if (anorm .ne. 0.0d0) rcond = ynorm/anorm
+ if (anorm .eq. 0.0d0) rcond = 0.0d0
+ 180 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dpodi.f b/com.oracle.truffle.r.native/src/dpodi.f
new file mode 100644
index 0000000000000000000000000000000000000000..4d52c5998fbefb4cf81542a8966d757e80d7650c
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dpodi.f
@@ -0,0 +1,122 @@
+c
+c dpodi computes the determinant and inverse of a certain
+c double precision symmetric positive definite matrix (see below)
+c using the factors computed by dpoco, dpofa or dqrdc.
+c
+c on entry
+c
+c a double precision(lda, n)
+c the output a from dpoco or dpofa
+c or the output x from dqrdc.
+c
+c lda integer
+c the leading dimension of the array a .
+c
+c n integer
+c the order of the matrix a .
+c
+c job integer
+c = 11 both determinant and inverse.
+c = 01 inverse only.
+c = 10 determinant only.
+c
+c on return
+c
+c a if dpoco or dpofa was used to factor a then
+c dpodi produces the upper half of inverse(a) .
+c if dqrdc was used to decompose x then
+c dpodi produces the upper half of inverse(trans(x)*x)
+c where trans(x) is the transpose.
+c elements of a below the diagonal are unchanged.
+c if the units digit of job is zero, a is unchanged.
+c
+c det double precision(2)
+c determinant of a or of trans(x)*x if requested.
+c otherwise not referenced.
+c determinant = det(1) * 10.0**det(2)
+c with 1.0 .le. det(1) .lt. 10.0
+c or det(1) .eq. 0.0 .
+c
+c error condition
+c
+c a division by zero will occur if the input factor contains
+c a zero on the diagonal and the inverse is requested.
+c it will not occur if the subroutines are called correctly
+c and if dpoco or dpofa has set info .eq. 0 .
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas daxpy,dscal
+c fortran mod
+c
+ subroutine dpodi(a,lda,n,det,job)
+ integer lda,n,job
+ double precision a(lda,*)
+ double precision det(2)
+c
+c internal variables
+c
+ double precision t
+ double precision s
+ integer i,j,jm1,k,kp1
+c
+c compute determinant
+c
+ if (job/10 .eq. 0) go to 70
+ det(1) = 1.0d0
+ det(2) = 0.0d0
+ s = 10.0d0
+ do 50 i = 1, n
+ det(1) = a(i,i)**2*det(1)
+c ...exit
+ if (det(1) .eq. 0.0d0) go to 60
+ 10 if (det(1) .ge. 1.0d0) go to 20
+ det(1) = s*det(1)
+ det(2) = det(2) - 1.0d0
+ go to 10
+ 20 continue
+ 30 if (det(1) .lt. s) go to 40
+ det(1) = det(1)/s
+ det(2) = det(2) + 1.0d0
+ go to 30
+ 40 continue
+ 50 continue
+ 60 continue
+ 70 continue
+c
+c compute inverse(r)
+c
+ if (mod(job,10) .eq. 0) go to 140
+ do 100 k = 1, n
+ a(k,k) = 1.0d0/a(k,k)
+ t = -a(k,k)
+ call dscal(k-1,t,a(1,k),1)
+ kp1 = k + 1
+ if (n .lt. kp1) go to 90
+ do 80 j = kp1, n
+ t = a(k,j)
+ a(k,j) = 0.0d0
+ call daxpy(k,t,a(1,k),1,a(1,j),1)
+ 80 continue
+ 90 continue
+ 100 continue
+c
+c form inverse(r) * trans(inverse(r))
+c
+ do 130 j = 1, n
+ jm1 = j - 1
+ if (jm1 .lt. 1) go to 120
+ do 110 k = 1, jm1
+ t = a(k,j)
+ call daxpy(k,t,a(1,j),1,a(1,k),1)
+ 110 continue
+ 120 continue
+ t = a(j,j)
+ call dscal(j,t,a(1,j),1)
+ 130 continue
+ 140 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dpofa.f b/com.oracle.truffle.r.native/src/dpofa.f
new file mode 100644
index 0000000000000000000000000000000000000000..94d71d20ede3f1788dd2bd6bf8ab9c11a23f6fef
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dpofa.f
@@ -0,0 +1,78 @@
+C Modified 2002-05-20 for R to add a tolerance of positive definiteness.
+C
+c
+c dpofa factors a double precision symmetric positive definite
+c matrix.
+c
+c dpofa is usually called by dpoco, but it can be called
+c directly with a saving in time if rcond is not needed.
+c (time for dpoco) = (1 + 18/n)*(time for dpofa) .
+c
+c on entry
+c
+c a double precision(lda, n)
+c the symmetric matrix to be factored. only the
+c diagonal and upper triangle are used.
+c
+c lda integer
+c the leading dimension of the array a .
+c
+c n integer
+c the order of the matrix a .
+c
+c on return
+c
+c a an upper triangular matrix r so that a = trans(r)*r
+c where trans(r) is the transpose.
+c the strict lower triangle is unaltered.
+c if info .ne. 0 , the factorization is not complete.
+c
+c info integer
+c = 0 for normal return.
+c = k signals an error condition. the leading minor
+c of order k is not positive definite.
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas ddot
+c fortran dsqrt
+c
+ subroutine dpofa(a,lda,n,info)
+ integer lda,n,info
+ double precision a(lda,*)
+c
+c internal variables
+c
+ double precision ddot,t,eps
+ double precision s
+ integer j,jm1,k
+ data eps/1.d-14/
+
+c begin block with ...exits to 40
+c
+c
+ do 30 j = 1, n
+ info = j
+ s = 0.0d0
+ jm1 = j - 1
+ if (jm1 .lt. 1) go to 20
+ do 10 k = 1, jm1
+ t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
+ t = t/a(k,k)
+ a(k,j) = t
+ s = s + t*t
+ 10 continue
+ 20 continue
+ s = a(j,j) - s
+c ......exit
+c if (s .le. 0.0d0) go to 40
+ if (s .le. eps * abs(a(j,j))) go to 40
+ a(j,j) = dsqrt(s)
+ 30 continue
+ info = 0
+ 40 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dposl.f b/com.oracle.truffle.r.native/src/dposl.f
new file mode 100644
index 0000000000000000000000000000000000000000..41c0e69502fc3cd94e3301364059ce3c1bdae1ac
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dposl.f
@@ -0,0 +1,72 @@
+c
+c dposl solves the double precision symmetric positive definite
+c system a * x = b
+c using the factors computed by dpoco or dpofa.
+c
+c on entry
+c
+c a double precision(lda, n)
+c the output from dpoco or dpofa.
+c
+c lda integer
+c the leading dimension of the array a .
+c
+c n integer
+c the order of the matrix a .
+c
+c b double precision(n)
+c the right hand side vector.
+c
+c on return
+c
+c b the solution vector x .
+c
+c error condition
+c
+c a division by zero will occur if the input factor contains
+c a zero on the diagonal. technically this indicates
+c singularity but it is usually caused by improper subroutine
+c arguments. it will not occur if the subroutines are called
+c correctly and info .eq. 0 .
+c
+c to compute inverse(a) * c where c is a matrix
+c with p columns
+c call dpoco(a,lda,n,rcond,z,info)
+c if (rcond is too small .or. info .ne. 0) go to ...
+c do 10 j = 1, p
+c call dposl(a,lda,n,c(1,j))
+c 10 continue
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas daxpy,ddot
+c
+ subroutine dposl(a,lda,n,b)
+ integer lda,n
+ double precision a(lda,*),b(*)
+c
+c internal variables
+c
+ double precision ddot,t
+ integer k,kb
+c
+c solve trans(r)*y = b
+c
+ do 10 k = 1, n
+ t = ddot(k-1,a(1,k),1,b(1),1)
+ b(k) = (b(k) - t)/a(k,k)
+ 10 continue
+c
+c solve r*x = y
+c
+ do 20 kb = 1, n
+ k = n + 1 - kb
+ b(k) = b(k)/a(k,k)
+ t = -b(k)
+ call daxpy(k-1,t,a(1,k),1,b(1),1)
+ 20 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dqrdc.f b/com.oracle.truffle.r.native/src/dqrdc.f
new file mode 100644
index 0000000000000000000000000000000000000000..54c9320cbea64ada95353a5231b39b158e168794
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dqrdc.f
@@ -0,0 +1,208 @@
+c
+c dqrdc uses householder transformations to compute the qr
+c factorization of an n by p matrix x. column pivoting
+c based on the 2-norms of the reduced columns may be
+c performed at the users option.
+c
+c on entry
+c
+c x double precision(ldx,p), where ldx .ge. n.
+c x contains the matrix whose decomposition is to be
+c computed.
+c
+c ldx integer.
+c ldx is the leading dimension of the array x.
+c
+c n integer.
+c n is the number of rows of the matrix x.
+c
+c p integer.
+c p is the number of columns of the matrix x.
+c
+c jpvt integer(p).
+c jpvt contains integers that control the selection
+c of the pivot columns. the k-th column x(k) of x
+c is placed in one of three classes according to the
+c value of jpvt(k).
+c
+c if jpvt(k) .gt. 0, then x(k) is an initial
+c column.
+c
+c if jpvt(k) .eq. 0, then x(k) is a free column.
+c
+c if jpvt(k) .lt. 0, then x(k) is a final column.
+c
+c before the decomposition is computed, initial columns
+c are moved to the beginning of the array x and final
+c columns to the end. both initial and final columns
+c are frozen in place during the computation and only
+c free columns are moved. at the k-th stage of the
+c reduction, if x(k) is occupied by a free column
+c it is interchanged with the free column of largest
+c reduced norm. jpvt is not referenced if
+c job .eq. 0.
+c
+c work double precision(p).
+c work is a work array. work is not referenced if
+c job .eq. 0.
+c
+c job integer.
+c job is an integer that initiates column pivoting.
+c if job .eq. 0, no pivoting is done.
+c if job .ne. 0, pivoting is done.
+c
+c on return
+c
+c x x contains in its upper triangle the upper
+c triangular matrix r of the qr factorization.
+c below its diagonal x contains information from
+c which the orthogonal part of the decomposition
+c can be recovered. note that if pivoting has
+c been requested, the decomposition is not that
+c of the original matrix x but that of x
+c with its columns permuted as described by jpvt.
+c
+c qraux double precision(p).
+c qraux contains further information required to recover
+c the orthogonal part of the decomposition.
+c
+c jpvt jpvt(k) contains the index of the column of the
+c original matrix that has been interchanged into
+c the k-th column, if pivoting was requested.
+c
+c linpack. this version dated 08/14/78 .
+c g.w. stewart, university of maryland, argonne national lab.
+c
+c dqrdc uses the following functions and subprograms.
+c
+c blas daxpy,ddot,dscal,dswap,dnrm2
+c fortran dabs,dmax1,min0,dsqrt
+c
+ subroutine dqrdc(x,ldx,n,p,qraux,jpvt,work,job)
+ integer ldx,n,p,job
+ integer jpvt(*)
+ double precision x(ldx,*),qraux(*),work(*)
+c
+c internal variables
+c
+ integer j,jp,jj, l,lp1,lup,maxj,pl,pu
+ double precision maxnrm,dnrm2,tt
+ double precision ddot,nrmxl,t
+ logical negj,swapj
+c
+c
+ pl = 1
+ pu = 0
+ if (job .eq. 0) go to 60
+c
+c pivoting has been requested. rearrange the columns
+c according to jpvt.
+c
+ do 20 j = 1, p
+ swapj = jpvt(j) .gt. 0
+ negj = jpvt(j) .lt. 0
+ jpvt(j) = j
+ if (negj) jpvt(j) = -j
+ if (.not.swapj) go to 10
+ if (j .ne. pl) call dswap(n,x(1,pl),1,x(1,j),1)
+ jpvt(j) = jpvt(pl)
+ jpvt(pl) = j
+ pl = pl + 1
+ 10 continue
+ 20 continue
+ pu = p
+ do 50 jj = 1, p
+ j = p - jj + 1
+ if (jpvt(j) .ge. 0) go to 40
+ jpvt(j) = -jpvt(j)
+ if (j .eq. pu) go to 30
+ call dswap(n,x(1,pu),1,x(1,j),1)
+ jp = jpvt(pu)
+ jpvt(pu) = jpvt(j)
+ jpvt(j) = jp
+ 30 continue
+ pu = pu - 1
+ 40 continue
+ 50 continue
+ 60 continue
+c
+c compute the norms of the free columns.
+c
+ if (pu .lt. pl) go to 80
+ do 70 j = pl, pu
+ qraux(j) = dnrm2(n,x(1,j),1)
+ work(j) = qraux(j)
+ 70 continue
+ 80 continue
+c
+c perform the householder reduction of x.
+c
+ lup = min0(n,p)
+ do 200 l = 1, lup
+ if (l .lt. pl .or. l .ge. pu) go to 120
+c
+c locate the column of largest norm and bring it
+c into the pivot position.
+c
+ maxnrm = 0.0d0
+ maxj = l
+ do 100 j = l, pu
+ if (qraux(j) .le. maxnrm) go to 90
+ maxnrm = qraux(j)
+ maxj = j
+ 90 continue
+ 100 continue
+ if (maxj .eq. l) go to 110
+ call dswap(n,x(1,l),1,x(1,maxj),1)
+ qraux(maxj) = qraux(l)
+ work(maxj) = work(l)
+ jp = jpvt(maxj)
+ jpvt(maxj) = jpvt(l)
+ jpvt(l) = jp
+ 110 continue
+ 120 continue
+ qraux(l) = 0.0d0
+ if (l .eq. n) go to 190
+c
+c compute the householder transformation for column l.
+c
+ nrmxl = dnrm2(n-l+1,x(l,l),1)
+ if (nrmxl .eq. 0.0d0) go to 180
+ if (x(l,l) .ne. 0.0d0) nrmxl = dsign(nrmxl,x(l,l))
+ call dscal(n-l+1,1.0d0/nrmxl,x(l,l),1)
+ x(l,l) = 1.0d0 + x(l,l)
+c
+c apply the transformation to the remaining columns,
+c updating the norms.
+c
+ lp1 = l + 1
+ if (p .lt. lp1) go to 170
+ do 160 j = lp1, p
+ t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
+ call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
+ if (j .lt. pl .or. j .gt. pu) go to 150
+ if (qraux(j) .eq. 0.0d0) go to 150
+ tt = 1.0d0 - (dabs(x(l,j))/qraux(j))**2
+ tt = dmax1(tt,0.0d0)
+ t = tt
+ tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j))**2
+ if (tt .eq. 1.0d0) go to 130
+ qraux(j) = qraux(j)*dsqrt(t)
+ go to 140
+ 130 continue
+ qraux(j) = dnrm2(n-l,x(l+1,j),1)
+ work(j) = qraux(j)
+ 140 continue
+ 150 continue
+ 160 continue
+ 170 continue
+c
+c save the transformation.
+c
+ qraux(l) = x(l,l)
+ x(l,l) = -nrmxl
+ 180 continue
+ 190 continue
+ 200 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dqrdc2.f b/com.oracle.truffle.r.native/src/dqrdc2.f
new file mode 100644
index 0000000000000000000000000000000000000000..6e9e23ff2d92bf847dbd2767ef3c426742b8a71a
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dqrdc2.f
@@ -0,0 +1,194 @@
+c
+c dqrdc2 uses householder transformations to compute the qr
+c factorization of an n by p matrix x. a limited column
+c pivoting strategy based on the 2-norms of the reduced columns
+c moves columns with near-zero norm to the right-hand edge of
+c the x matrix. this strategy means that sequential one
+c degree-of-freedom effects can be computed in a natural way.
+c
+c i am very nervous about modifying linpack code in this way.
+c if you are a computational linear algebra guru and you really
+c understand how to solve this problem please feel free to
+c suggest improvements to this code.
+c
+c Another change was to compute the rank.
+c
+c on entry
+c
+c x double precision(ldx,p), where ldx .ge. n.
+c x contains the matrix whose decomposition is to be
+c computed.
+c
+c ldx integer.
+c ldx is the leading dimension of the array x.
+c
+c n integer.
+c n is the number of rows of the matrix x.
+c
+c p integer.
+c p is the number of columns of the matrix x.
+c
+c tol double precision
+c tol is the nonnegative tolerance used to
+c determine the subset of the columns of x
+c included in the solution.
+c
+c jpvt integer(p).
+c integers which are swapped in the same way as the
+c the columns of x during pivoting. on entry these
+c should be set equal to the column indices of the
+c columns of the x matrix (typically 1 to p).
+c
+c work double precision(p,2).
+c work is a work array.
+c
+c on return
+c
+c x x contains in its upper triangle the upper
+c triangular matrix r of the qr factorization.
+c below its diagonal x contains information from
+c which the orthogonal part of the decomposition
+c can be recovered. note that if pivoting has
+c been requested, the decomposition is not that
+c of the original matrix x but that of x
+c with its columns permuted as described by jpvt.
+c
+c k integer.
+c k contains the number of columns of x judged
+c to be linearly independent.
+c
+c qraux double precision(p).
+c qraux contains further information required to recover
+c the orthogonal part of the decomposition.
+c
+c jpvt jpvt(k) contains the index of the column of the
+c original matrix that has been interchanged into
+c the k-th column.
+c
+c original (dqrdc.f) linpack version dated 08/14/78 .
+c g.w. stewart, university of maryland, argonne national lab.
+c
+c this version dated 22 august 1995
+c ross ihaka
+c
+c bug fixes 29 September 1999 BDR (p > n case, inaccurate ranks)
+c
+c
+c dqrdc2 uses the following functions and subprograms.
+c
+c blas daxpy,ddot,dscal,dnrm2
+c fortran dabs,dmax1,min0,dsqrt
+c
+ subroutine dqrdc2(x,ldx,n,p,tol,k,qraux,jpvt,work)
+ integer ldx,n,p
+ integer jpvt(*)
+ double precision x(ldx,*),qraux(*),work(p,2),tol
+c
+c internal variables
+c
+ integer i,j,l,lp1,lup,k
+ double precision dnrm2,tt,ttt
+ double precision ddot,nrmxl,t
+c
+c
+c compute the norms of the columns of x.
+c
+ do 70 j = 1, p
+ qraux(j) = dnrm2(n,x(1,j),1)
+ work(j,1) = qraux(j)
+ work(j,2) = qraux(j)
+ if(work(j,2) .eq. 0.0d0) work(j,2) = 1.0d0
+ 70 continue
+c
+c perform the householder reduction of x.
+c
+ lup = min0(n,p)
+ k = p + 1
+ do 200 l = 1, lup
+c
+c previous version only cycled l to lup
+c
+c cycle the columns from l to p left-to-right until one
+c with non-negligible norm is located. a column is considered
+c to have become negligible if its norm has fallen below
+c tol times its original norm. the check for l .le. k
+c avoids infinite cycling.
+c
+ 80 continue
+ if (l .ge. k .or. qraux(l) .ge. work(l,2)*tol) go to 120
+ lp1 = l+1
+ do 100 i=1,n
+ t = x(i,l)
+ do 90 j=lp1,p
+ x(i,j-1) = x(i,j)
+ 90 continue
+ x(i,p) = t
+ 100 continue
+ i = jpvt(l)
+ t = qraux(l)
+ tt = work(l,1)
+ ttt = work(l,2)
+ do 110 j=lp1,p
+ jpvt(j-1) = jpvt(j)
+ qraux(j-1) = qraux(j)
+ work(j-1,1) = work(j,1)
+ work(j-1,2) = work(j,2)
+ 110 continue
+ jpvt(p) = i
+ qraux(p) = t
+ work(p,1) = tt
+ work(p,2) = ttt
+ k = k - 1
+ go to 80
+ 120 continue
+ if (l .eq. n) go to 190
+c
+c compute the householder transformation for column l.
+c
+ nrmxl = dnrm2(n-l+1,x(l,l),1)
+ if (nrmxl .eq. 0.0d0) go to 180
+ if (x(l,l) .ne. 0.0d0) nrmxl = dsign(nrmxl,x(l,l))
+ call dscal(n-l+1,1.0d0/nrmxl,x(l,l),1)
+ x(l,l) = 1.0d0 + x(l,l)
+c
+c apply the transformation to the remaining columns,
+c updating the norms.
+c
+ lp1 = l + 1
+ if (p .lt. lp1) go to 170
+ do 160 j = lp1, p
+ t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
+ call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
+ if (qraux(j) .eq. 0.0d0) go to 150
+ tt = 1.0d0 - (dabs(x(l,j))/qraux(j))**2
+ tt = dmax1(tt,0.0d0)
+ t = tt
+c
+c modified 9/99 by BDR. Re-compute norms if there is large reduction
+c The tolerance here is on the squared norm
+c In this version we need accurate norms, so re-compute often.
+c work(j,1) is only updated in one case: looks like a bug -- no longer used
+c
+c tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j,1))**2
+c if (tt .eq. 1.0d0) go to 130
+ if (dabs(t) .lt. 1d-6) go to 130
+ qraux(j) = qraux(j)*dsqrt(t)
+ go to 140
+ 130 continue
+ qraux(j) = dnrm2(n-l,x(l+1,j),1)
+ work(j,1) = qraux(j)
+ 140 continue
+ 150 continue
+ 160 continue
+ 170 continue
+c
+c save the transformation.
+c
+ qraux(l) = x(l,l)
+ x(l,l) = -nrmxl
+ 180 continue
+ 190 continue
+ 200 continue
+ k = min0(k - 1, n)
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dqrls.f b/com.oracle.truffle.r.native/src/dqrls.f
new file mode 100644
index 0000000000000000000000000000000000000000..31e029f5a02fe6201476b3c26b24fb9a915c67f7
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dqrls.f
@@ -0,0 +1,116 @@
+c
+c dqrfit is a subroutine to compute least squares solutions
+c to the system
+c
+c (1) x * b = y
+c
+c which may be either under-determined or over-determined.
+c the user must supply a tolerance to limit the columns of
+c x used in computing the solution. in effect, a set of
+c columns with a condition number approximately bounded by
+c 1/tol is used, the other components of b being set to zero.
+c
+c on entry
+c
+c x double precision(n,p).
+c x contains n-by-p coefficient matrix of
+c the system (1), x is destroyed by dqrfit.
+c
+c n the number of rows of the matrix x.
+c
+c p the number of columns of the matrix x.
+c
+c y double precision(n,ny)
+c y contains the right hand side(s) of the system (1).
+c
+c ny the number of right hand sides of the system (1).
+c
+c tol double precision
+c tol is the nonnegative tolerance used to
+c determine the subset of columns of x included
+c in the solution. columns are pivoted out of
+c decomposition if
+c
+c jpvt integer(p)
+c the values in jpvt are permuted in the same
+c way as the columns of x. this can be useful
+c in unscrambling coefficients etc.
+c
+c work double precision(2*p)
+c work is an array used by dqrdc2 and dqrsl.
+c
+c on return
+c
+c x contains the output array from dqrdc2.
+c namely the qr decomposition of x stored in
+c compact form.
+c
+c b double precision(p,ny)
+c b contains the solution vectors with rows permuted
+c in the same way as the columns of x. components
+c corresponding to columns not used are set to zero.
+c
+c rsd double precision(n,ny)
+c rsd contains the residual vectors y-x*b.
+c
+c qty double precision(n,ny) t
+c qty contains the vectors q y. note that
+c the initial p elements of this vector are
+c permuted in the same way as the columns of x.
+c
+c k integer
+c k contains the number of columns used in the
+c solution.
+c
+c jpvt has its contents permuted as described above.
+c
+c qraux double precision(p)
+c qraux contains auxiliary information on the
+c qr decomposition of x.
+c
+c
+c on return the arrays x, jpvt and qraux contain the
+c usual output from dqrdc, so that the qr decomposition
+c of x with pivoting is fully available to the user.
+c in particular, columns jpvt(1), jpvt(2),...,jpvt(k)
+c were used in the solution, and the condition number
+c associated with those columns is estimated by
+c abs(x(1,1)/x(k,k)).
+c
+c dqrfit uses the linpack routines dqrdc and dqrsl.
+c
+ subroutine dqrls(x,n,p,y,ny,tol,b,rsd,qty,k,jpvt,qraux,work)
+ integer n,p,ny,k,jpvt(p)
+ double precision x(n,p),y(n,ny),tol,b(p,ny),rsd(n,ny),
+ . qty(n,ny),qraux(p),work(p)
+c
+c internal variables.
+c
+ integer info,j,jj,kk
+c
+c reduce x.
+c
+ call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)
+c
+c solve the truncated least squares problem for each rhs.
+c
+ if(k .gt. 0) then
+ do 20 jj=1,ny
+ 20 call dqrsl(x,n,n,k,qraux,y(1,jj),rsd(1,jj),qty(1,jj),
+ 1 b(1,jj),rsd(1,jj),rsd(1,jj),1110,info)
+ else
+ do 30 i=1,n
+ do 30 jj=1,ny
+ 30 rsd(i,jj) = y(i,jj)
+ endif
+c
+c set the unused components of b to zero.
+c
+ kk = k + 1
+ do 50 j=kk,p
+ do 40 jj=1,ny
+ b(j,jj) = 0.d0
+ 40 continue
+ 50 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dqrsl.f b/com.oracle.truffle.r.native/src/dqrsl.f
new file mode 100644
index 0000000000000000000000000000000000000000..aab138a5feb558905e7f4c597afdb2c1b14002a8
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dqrsl.f
@@ -0,0 +1,275 @@
+c
+c dqrsl applies the output of dqrdc to compute coordinate
+c transformations, projections, and least squares solutions.
+c for k .le. min(n,p), let xk be the matrix
+c
+c xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k)))
+c
+c formed from columnns jpvt(1), ... ,jpvt(k) of the original
+c n x p matrix x that was input to dqrdc (if no pivoting was
+c done, xk consists of the first k columns of x in their
+c original order). dqrdc produces a factored orthogonal matrix q
+c and an upper triangular matrix r such that
+c
+c xk = q * (r)
+c (0)
+c
+c this information is contained in coded form in the arrays
+c x and qraux.
+c
+c on entry
+c
+c x double precision(ldx,p).
+c x contains the output of dqrdc.
+c
+c ldx integer.
+c ldx is the leading dimension of the array x.
+c
+c n integer.
+c n is the number of rows of the matrix xk. it must
+c have the same value as n in dqrdc.
+c
+c k integer.
+c k is the number of columns of the matrix xk. k
+c must nnot be greater than min(n,p), where p is the
+c same as in the calling sequence to dqrdc.
+c
+c qraux double precision(p).
+c qraux contains the auxiliary output from dqrdc.
+c
+c y double precision(n)
+c y contains an n-vector that is to be manipulated
+c by dqrsl.
+c
+c job integer.
+c job specifies what is to be computed. job has
+c the decimal expansion abcde, with the following
+c meaning.
+c
+c if a.ne.0, compute qy.
+c if b,c,d, or e .ne. 0, compute qty.
+c if c.ne.0, compute b.
+c if d.ne.0, compute rsd.
+c if e.ne.0, compute xb.
+c
+c note that a request to compute b, rsd, or xb
+c automatically triggers the computation of qty, for
+c which an array must be provided in the calling
+c sequence.
+c
+c on return
+c
+c qy double precision(n).
+c qy conntains q*y, if its computation has been
+c requested.
+c
+c qty double precision(n).
+c qty contains trans(q)*y, if its computation has
+c been requested. here trans(q) is the
+c transpose of the matrix q.
+c
+c b double precision(k)
+c b contains the solution of the least squares problem
+c
+c minimize norm2(y - xk*b),
+c
+c if its computation has been requested. (note that
+c if pivoting was requested in dqrdc, the j-th
+c component of b will be associated with column jpvt(j)
+c of the original matrix x that was input into dqrdc.)
+c
+c rsd double precision(n).
+c rsd contains the least squares residual y - xk*b,
+c if its computation has been requested. rsd is
+c also the orthogonal projection of y onto the
+c orthogonal complement of the column space of xk.
+c
+c xb double precision(n).
+c xb contains the least squares approximation xk*b,
+c if its computation has been requested. xb is also
+c the orthogonal projection of y onto the column space
+c of x.
+c
+c info integer.
+c info is zero unless the computation of b has
+c been requested and r is exactly singular. in
+c this case, info is the index of the first zero
+c diagonal element of r and b is left unaltered.
+c
+c the parameters qy, qty, b, rsd, and xb are not referenced
+c if their computation is not requested and in this case
+c can be replaced by dummy variables in the calling program.
+c to save storage, the user may in some cases use the same
+c array for different parameters in the calling sequence. a
+c frequently occuring example is when one wishes to compute
+c any of b, rsd, or xb and does not need y or qty. in this
+c case one may identify y, qty, and one of b, rsd, or xb, while
+c providing separate arrays for anything else that is to be
+c computed. thus the calling sequence
+c
+c call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info)
+c
+c will result in the computation of b and rsd, with rsd
+c overwriting y. more generally, each item in the following
+c list contains groups of permissible identifications for
+c a single callinng sequence.
+c
+c 1. (y,qty,b) (rsd) (xb) (qy)
+c
+c 2. (y,qty,rsd) (b) (xb) (qy)
+c
+c 3. (y,qty,xb) (b) (rsd) (qy)
+c
+c 4. (y,qy) (qty,b) (rsd) (xb)
+c
+c 5. (y,qy) (qty,rsd) (b) (xb)
+c
+c 6. (y,qy) (qty,xb) (b) (rsd)
+c
+c in any group the value returned in the array allocated to
+c the group corresponds to the last member of the group.
+c
+c linpack. this version dated 08/14/78 .
+c g.w. stewart, university of maryland, argonne national lab.
+c
+c dqrsl uses the following functions and subprograms.
+c
+c BLAS daxpy,dcopy,ddot
+c Fortran dabs,min0,mod
+c
+ subroutine dqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
+ integer ldx,n,k,job,info
+ double precision x(ldx,*),qraux(*),y(*),qy(*),qty(*),b(*),rsd(*),
+ * xb(*)
+c
+c internal variables
+c
+ integer i,j,jj,ju,kp1
+ double precision ddot,t,temp
+ logical cb,cqy,cqty,cr,cxb
+c
+c
+c set info flag.
+c
+ info = 0
+c
+c determine what is to be computed.
+c
+ cqy = job/10000 .ne. 0
+ cqty = mod(job,10000) .ne. 0
+ cb = mod(job,1000)/100 .ne. 0
+ cr = mod(job,100)/10 .ne. 0
+ cxb = mod(job,10) .ne. 0
+ ju = min0(k,n-1)
+c
+c special action when n=1.
+c
+ if (ju .ne. 0) go to 40
+ if (cqy) qy(1) = y(1)
+ if (cqty) qty(1) = y(1)
+ if (cxb) xb(1) = y(1)
+ if (.not.cb) go to 30
+ if (x(1,1) .ne. 0.0d0) go to 10
+ info = 1
+ go to 20
+ 10 continue
+ b(1) = y(1)/x(1,1)
+ 20 continue
+ 30 continue
+ if (cr) rsd(1) = 0.0d0
+ go to 250
+ 40 continue
+c
+c set up to compute qy or qty.
+c
+ if (cqy) call dcopy(n,y,1,qy,1)
+ if (cqty) call dcopy(n,y,1,qty,1)
+ if (.not.cqy) go to 70
+c
+c compute qy.
+c
+ do 60 jj = 1, ju
+ j = ju - jj + 1
+ if (qraux(j) .eq. 0.0d0) go to 50
+ temp = x(j,j)
+ x(j,j) = qraux(j)
+ t = -ddot(n-j+1,x(j,j),1,qy(j),1)/x(j,j)
+ call daxpy(n-j+1,t,x(j,j),1,qy(j),1)
+ x(j,j) = temp
+ 50 continue
+ 60 continue
+ 70 continue
+ if (.not.cqty) go to 100
+c
+c compute trans(q)*y.
+c
+ do 90 j = 1, ju
+ if (qraux(j) .eq. 0.0d0) go to 80
+ temp = x(j,j)
+ x(j,j) = qraux(j)
+ t = -ddot(n-j+1,x(j,j),1,qty(j),1)/x(j,j)
+ call daxpy(n-j+1,t,x(j,j),1,qty(j),1)
+ x(j,j) = temp
+ 80 continue
+ 90 continue
+ 100 continue
+c
+c set up to compute b, rsd, or xb.
+c
+ if (cb) call dcopy(k,qty,1,b,1)
+ kp1 = k + 1
+ if (cxb) call dcopy(k,qty,1,xb,1)
+ if (cr .and. k .lt. n) call dcopy(n-k,qty(kp1),1,rsd(kp1),1)
+ if (.not.cxb .or. kp1 .gt. n) go to 120
+ do 110 i = kp1, n
+ xb(i) = 0.0d0
+ 110 continue
+ 120 continue
+ if (.not.cr) go to 140
+ do 130 i = 1, k
+ rsd(i) = 0.0d0
+ 130 continue
+ 140 continue
+ if (.not.cb) go to 190
+c
+c compute b.
+c
+ do 170 jj = 1, k
+ j = k - jj + 1
+ if (x(j,j) .ne. 0.0d0) go to 150
+ info = j
+c ......exit
+ go to 180
+ 150 continue
+ b(j) = b(j)/x(j,j)
+ if (j .eq. 1) go to 160
+ t = -b(j)
+ call daxpy(j-1,t,x(1,j),1,b,1)
+ 160 continue
+ 170 continue
+ 180 continue
+ 190 continue
+ if (.not.cr .and. .not.cxb) go to 240
+c
+c compute rsd or xb as required.
+c
+ do 230 jj = 1, ju
+ j = ju - jj + 1
+ if (qraux(j) .eq. 0.0d0) go to 220
+ temp = x(j,j)
+ x(j,j) = qraux(j)
+ if (.not.cr) go to 200
+ t = -ddot(n-j+1,x(j,j),1,rsd(j),1)/x(j,j)
+ call daxpy(n-j+1,t,x(j,j),1,rsd(j),1)
+ 200 continue
+ if (.not.cxb) go to 210
+ t = -ddot(n-j+1,x(j,j),1,xb(j),1)/x(j,j)
+ call daxpy(n-j+1,t,x(j,j),1,xb(j),1)
+ 210 continue
+ x(j,j) = temp
+ 220 continue
+ 230 continue
+ 240 continue
+ 250 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dqrutl.f b/com.oracle.truffle.r.native/src/dqrutl.f
new file mode 100644
index 0000000000000000000000000000000000000000..2d208ca0d8cb33e6215801c7637bc9465f167181
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dqrutl.f
@@ -0,0 +1,66 @@
+c dqr Utilities: Interface to the different "switches" of dqrsl().
+c
+ subroutine dqrqty(x, n, k, qraux, y, ny, qty)
+
+ integer n, k, ny
+ double precision x(n,k), qraux(k), y(n,ny), qty(n,ny)
+ integer info, j
+ double precision dummy(1)
+ do 10 j = 1,ny
+ call dqrsl(x, n, n, k, qraux, y(1,j), dummy, qty(1,j),
+ & dummy, dummy, dummy, 1000, info)
+ 10 continue
+ return
+ end
+c
+ subroutine dqrqy(x, n, k, qraux, y, ny, qy)
+
+ integer n, k, ny
+ double precision x(n,k), qraux(k), y(n,ny), qy(n,ny)
+ integer info, j
+ double precision dummy(1)
+ do 10 j = 1,ny
+ call dqrsl(x, n, n, k, qraux, y(1,j), qy(1,j),
+ & dummy, dummy, dummy, dummy, 10000, info)
+ 10 continue
+ return
+ end
+c
+ subroutine dqrcf(x, n, k, qraux, y, ny, b, info)
+
+ integer n, k, ny, info
+ double precision x(n,k), qraux(k), y(n,ny), b(k,ny)
+ integer j
+ double precision dummy(1)
+ do 10 j = 1,ny
+ call dqrsl(x, n, n, k, qraux, y(1,j), dummy,
+ & y(1,j), b(1,j), dummy, dummy, 100, info)
+ 10 continue
+ return
+ end
+c
+ subroutine dqrrsd(x, n, k, qraux, y, ny, rsd)
+
+ integer n, k, ny
+ double precision x(n,k), qraux(k), y(n,ny), rsd(n,ny)
+ integer info, j
+ double precision dummy(1)
+ do 10 j = 1,ny
+ call dqrsl(x, n, n, k, qraux, y(1,j), dummy,
+ & y(1,j), dummy, rsd(1,j), dummy, 10, info)
+ 10 continue
+ return
+ end
+c
+ subroutine dqrxb(x, n, k, qraux, y, ny, xb)
+
+ integer n, k, ny
+ double precision x(n,k), qraux(k), y(n,ny), xb(n,ny)
+ integer info, j
+ double precision dummy(1)
+ do 10 j = 1,ny
+ call dqrsl(x, n, n, k, qraux, y(1,j), dummy,
+ & y(1,j), dummy, dummy, xb(1,j), 1, info)
+ 10 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dsvdc.f b/com.oracle.truffle.r.native/src/dsvdc.f
new file mode 100644
index 0000000000000000000000000000000000000000..bdfd97d77a42a9a378a9e2b06ab802e325ee9505
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dsvdc.f
@@ -0,0 +1,495 @@
+c -- called from R's svd(x, ..., LINPACK = TRUE) , i.e, *NOT* by default --
+c
+c dsvdc is a subroutine to reduce a double precision nxp matrix x
+c by orthogonal transformations u and v to diagonal form. the
+c diagonal elements s(i) are the singular values of x. the
+c columns of u are the corresponding left singular vectors,
+c and the columns of v the right singular vectors.
+c
+c on entry
+c
+c x double precision(ldx,p), where ldx.ge.n.
+c x contains the matrix whose singular value
+c decomposition is to be computed. x is
+c destroyed by dsvdc.
+c
+c ldx integer.
+c ldx is the leading dimension of the array x.
+c
+c n integer.
+c n is the number of rows of the matrix x.
+c
+c p integer.
+c p is the number of columns of the matrix x.
+c
+c ldu integer.
+c ldu is the leading dimension of the array u.
+c (see below).
+c
+c ldv integer.
+c ldv is the leading dimension of the array v.
+c (see below).
+c
+c work double precision(n).
+c work is a scratch array.
+c
+c job integer.
+c job controls the computation of the singular
+c vectors. it has the decimal expansion ab
+c with the following meaning
+c
+c a.eq.0 do not compute the left singular
+c vectors.
+c a.eq.1 return the n left singular vectors
+c in u.
+c a.ge.2 return the first min(n,p) singular
+c vectors in u.
+c b.eq.0 do not compute the right singular
+c vectors.
+c b.eq.1 return the right singular vectors
+c in v.
+c
+c on return
+c
+c s double precision(mm), where mm=min(n+1,p).
+c the first min(n,p) entries of s contain the
+c singular values of x arranged in descending
+c order of magnitude.
+c
+c e double precision(p),
+c e ordinarily contains zeros. however see the
+c discussion of info for exceptions.
+c
+c u double precision(ldu,k), where ldu.ge.n. if
+c joba.eq.1 then k.eq.n, if joba.ge.2
+c then k.eq.min(n,p).
+c u contains the matrix of left singular vectors.
+c u is not referenced if joba.eq.0. if n.le.p
+c or if joba.eq.2, then u may be identified with x
+c in the subroutine call.
+c
+c v double precision(ldv,p), where ldv.ge.p.
+c v contains the matrix of right singular vectors.
+c v is not referenced if job.eq.0. if p.le.n,
+c then v may be identified with x in the
+c subroutine call.
+c
+c info integer.
+c the singular values (and their corresponding
+c singular vectors) s(info+1),s(info+2),...,s(m)
+c are correct (here m=min(n,p)). thus if
+c info.eq.0, all the singular values and their
+c vectors are correct. in any event, the matrix
+c b = trans(u)*x*v is the bidiagonal matrix
+c with the elements of s on its diagonal and the
+c elements of e on its super-diagonal (trans(u)
+c is the transpose of u). thus the singular
+c values of x and b are the same.
+c
+c linpack. this version dated 08/14/78 .
+c correction made to shift 2/84.
+c g.w. stewart, university of maryland, argonne national lab.
+c
+c Modified 2000-12-28 to use a relative convergence test,
+c as this was infinite-looping on ix86.
+c
+c dsvdc uses the following functions and subprograms.
+c
+c external drot
+c blas daxpy,ddot,dscal,dswap,dnrm2,drotg
+c fortran dabs,dmax1,max0,min0,mod,dsqrt
+c
+ subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
+ integer ldx,n,p,ldu,ldv,job,info
+ double precision x(ldx,*),s(*),e(*),u(ldu,*),v(ldv,*),work(*)
+c
+c internal variables
+c
+ integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit,
+ * mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1
+ double precision ddot,t
+ double precision b,c,cs,el,emm1,f,g,dnrm2,scale,shift,sl,sm,sn,
+ * smm1,t1,test,ztest,acc
+ logical wantu,wantv
+c
+c unnecessary initializations of l and ls to keep g77 -Wall happy
+c
+ l = 0
+ ls = 0
+c
+c
+c set the maximum number of iterations.
+c
+ maxit = 30
+c
+c determine what is to be computed.
+c
+ wantu = .false.
+ wantv = .false.
+ jobu = mod(job,100)/10
+ ncu = n
+ if (jobu .gt. 1) ncu = min0(n,p)
+ if (jobu .ne. 0) wantu = .true.
+ if (mod(job,10) .ne. 0) wantv = .true.
+c
+c reduce x to bidiagonal form, storing the diagonal elements
+c in s and the super-diagonal elements in e.
+c
+ info = 0
+ nct = min0(n-1,p)
+ nrt = max0(0,min0(p-2,n))
+ lu = max0(nct,nrt)
+ if (lu .lt. 1) go to 170
+ do 160 l = 1, lu
+ lp1 = l + 1
+ if (l .gt. nct) go to 20
+c
+c compute the transformation for the l-th column and
+c place the l-th diagonal in s(l).
+c
+ s(l) = dnrm2(n-l+1,x(l,l),1)
+ if (s(l) .eq. 0.0d0) go to 10
+ if (x(l,l) .ne. 0.0d0) s(l) = dsign(s(l),x(l,l))
+ call dscal(n-l+1,1.0d0/s(l),x(l,l),1)
+ x(l,l) = 1.0d0 + x(l,l)
+ 10 continue
+ s(l) = -s(l)
+ 20 continue
+ if (p .lt. lp1) go to 50
+ do 40 j = lp1, p
+ if (l .gt. nct) go to 30
+ if (s(l) .eq. 0.0d0) go to 30
+c
+c apply the transformation.
+c
+ t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
+ call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
+ 30 continue
+c
+c place the l-th row of x into e for the
+c subsequent calculation of the row transformation.
+c
+ e(j) = x(l,j)
+ 40 continue
+ 50 continue
+ if (.not.wantu .or. l .gt. nct) go to 70
+c
+c place the transformation in u for subsequent back
+c multiplication.
+c
+ do 60 i = l, n
+ u(i,l) = x(i,l)
+ 60 continue
+ 70 continue
+ if (l .gt. nrt) go to 150
+c
+c compute the l-th row transformation and place the
+c l-th super-diagonal in e(l).
+c
+ e(l) = dnrm2(p-l,e(lp1),1)
+ if (e(l) .eq. 0.0d0) go to 80
+ if (e(lp1) .ne. 0.0d0) e(l) = dsign(e(l),e(lp1))
+ call dscal(p-l,1.0d0/e(l),e(lp1),1)
+ e(lp1) = 1.0d0 + e(lp1)
+ 80 continue
+ e(l) = -e(l)
+ if (lp1 .gt. n .or. e(l) .eq. 0.0d0) go to 120
+c
+c apply the transformation.
+c
+ do 90 i = lp1, n
+ work(i) = 0.0d0
+ 90 continue
+ do 100 j = lp1, p
+ call daxpy(n-l,e(j),x(lp1,j),1,work(lp1),1)
+ 100 continue
+ do 110 j = lp1, p
+ call daxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1)
+ 110 continue
+ 120 continue
+ if (.not.wantv) go to 140
+c
+c place the transformation in v for subsequent
+c back multiplication.
+c
+ do 130 i = lp1, p
+ v(i,l) = e(i)
+ 130 continue
+ 140 continue
+ 150 continue
+ 160 continue
+ 170 continue
+c
+c set up the final bidiagonal matrix or order m.
+c
+ m = min0(p,n+1)
+ nctp1 = nct + 1
+ nrtp1 = nrt + 1
+ if (nct .lt. p) s(nctp1) = x(nctp1,nctp1)
+ if (n .lt. m) s(m) = 0.0d0
+ if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m)
+ e(m) = 0.0d0
+c
+c if required, generate u.
+c
+ if (.not.wantu) go to 300
+ if (ncu .lt. nctp1) go to 200
+ do 190 j = nctp1, ncu
+ do 180 i = 1, n
+ u(i,j) = 0.0d0
+ 180 continue
+ u(j,j) = 1.0d0
+ 190 continue
+ 200 continue
+ if (nct .lt. 1) go to 290
+ do 280 ll = 1, nct
+ l = nct - ll + 1
+ if (s(l) .eq. 0.0d0) go to 250
+ lp1 = l + 1
+ if (ncu .lt. lp1) go to 220
+ do 210 j = lp1, ncu
+ t = -ddot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l)
+ call daxpy(n-l+1,t,u(l,l),1,u(l,j),1)
+ 210 continue
+ 220 continue
+ call dscal(n-l+1,-1.0d0,u(l,l),1)
+ u(l,l) = 1.0d0 + u(l,l)
+ lm1 = l - 1
+ if (lm1 .lt. 1) go to 240
+ do 230 i = 1, lm1
+ u(i,l) = 0.0d0
+ 230 continue
+ 240 continue
+ go to 270
+ 250 continue
+ do 260 i = 1, n
+ u(i,l) = 0.0d0
+ 260 continue
+ u(l,l) = 1.0d0
+ 270 continue
+ 280 continue
+ 290 continue
+ 300 continue
+c
+c if it is required, generate v.
+c
+ if (.not.wantv) go to 350
+ do 340 ll = 1, p
+ l = p - ll + 1
+ lp1 = l + 1
+ if (l .gt. nrt) go to 320
+ if (e(l) .eq. 0.0d0) go to 320
+ do 310 j = lp1, p
+ t = -ddot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l)
+ call daxpy(p-l,t,v(lp1,l),1,v(lp1,j),1)
+ 310 continue
+ 320 continue
+ do 330 i = 1, p
+ v(i,l) = 0.0d0
+ 330 continue
+ v(l,l) = 1.0d0
+ 340 continue
+ 350 continue
+c
+c main iteration loop for the singular values.
+c
+ mm = m
+ iter = 0
+ 360 continue
+c
+c quit if all the singular values have been found.
+c
+c ...exit
+ if (m .eq. 0) go to 620
+c
+c if too many iterations have been performed, set
+c flag and return.
+c
+ if (iter .lt. maxit) go to 370
+ info = m
+c ......exit
+ go to 620
+ 370 continue
+c
+c this section of the program inspects for
+c negligible elements in the s and e arrays. on
+c completion the variables kase and l are set as follows.
+c
+c kase = 1 if s(m) and e(l-1) are negligible and l.lt.m
+c kase = 2 if s(l) is negligible and l.lt.m
+c kase = 3 if e(l-1) is negligible, l.lt.m, and
+c s(l), ..., s(m) are not negligible (qr step).
+c kase = 4 if e(m-1) is negligible (convergence).
+c
+ do 390 ll = 1, m
+ l = m - ll
+c ...exit
+ if (l .eq. 0) go to 400
+ test = dabs(s(l)) + dabs(s(l+1))
+ ztest = test + dabs(e(l))
+ acc = dabs(test - ztest)/(1.0d-100 + test)
+ if (acc .gt. 1.d-15) goto 380
+c if (ztest .ne. test) go to 380
+ e(l) = 0.0d0
+c ......exit
+ go to 400
+ 380 continue
+ 390 continue
+ 400 continue
+ if (l .ne. m - 1) go to 410
+ kase = 4
+ go to 480
+ 410 continue
+ lp1 = l + 1
+ mp1 = m + 1
+ do 430 lls = lp1, mp1
+ ls = m - lls + lp1
+c ...exit
+ if (ls .eq. l) go to 440
+ test = 0.0d0
+ if (ls .ne. m) test = test + dabs(e(ls))
+ if (ls .ne. l + 1) test = test + dabs(e(ls-1))
+ ztest = test + dabs(s(ls))
+c 1.0d-100 is to guard against a zero matrix, hence zero test
+ acc = dabs(test - ztest)/(1.0d-100 + test)
+ if (acc .gt. 1.d-15) goto 420
+c if (ztest .ne. test) go to 420
+ s(ls) = 0.0d0
+c ......exit
+ go to 440
+ 420 continue
+ 430 continue
+ 440 continue
+ if (ls .ne. l) go to 450
+ kase = 3
+ go to 470
+ 450 continue
+ if (ls .ne. m) go to 460
+ kase = 1
+ go to 470
+ 460 continue
+ kase = 2
+ l = ls
+ 470 continue
+ 480 continue
+ l = l + 1
+c
+c perform the task indicated by kase.
+c
+ go to (490,520,540,570), kase
+c
+c deflate negligible s(m).
+c
+ 490 continue
+ mm1 = m - 1
+ f = e(m-1)
+ e(m-1) = 0.0d0
+ do 510 kk = l, mm1
+ k = mm1 - kk + l
+ t1 = s(k)
+ call drotg(t1,f,cs,sn)
+ s(k) = t1
+ if (k .eq. l) go to 500
+ f = -sn*e(k-1)
+ e(k-1) = cs*e(k-1)
+ 500 continue
+ if (wantv) call drot(p,v(1,k),1,v(1,m),1,cs,sn)
+ 510 continue
+ go to 610
+c
+c split at negligible s(l).
+c
+ 520 continue
+ f = e(l-1)
+ e(l-1) = 0.0d0
+ do 530 k = l, m
+ t1 = s(k)
+ call drotg(t1,f,cs,sn)
+ s(k) = t1
+ f = -sn*e(k)
+ e(k) = cs*e(k)
+ if (wantu) call drot(n,u(1,k),1,u(1,l-1),1,cs,sn)
+ 530 continue
+ go to 610
+c
+c perform one qr step.
+c
+ 540 continue
+c
+c calculate the shift.
+c
+ scale = dmax1(dabs(s(m)),dabs(s(m-1)),dabs(e(m-1)),
+ * dabs(s(l)),dabs(e(l)))
+ sm = s(m)/scale
+ smm1 = s(m-1)/scale
+ emm1 = e(m-1)/scale
+ sl = s(l)/scale
+ el = e(l)/scale
+ b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0
+ c = (sm*emm1)**2
+ shift = 0.0d0
+ if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550
+ shift = dsqrt(b**2+c)
+ if (b .lt. 0.0d0) shift = -shift
+ shift = c/(b + shift)
+ 550 continue
+ f = (sl + sm)*(sl - sm) + shift
+ g = sl*el
+c
+c chase zeros.
+c
+ mm1 = m - 1
+ do 560 k = l, mm1
+ call drotg(f,g,cs,sn)
+ if (k .ne. l) e(k-1) = f
+ f = cs*s(k) + sn*e(k)
+ e(k) = cs*e(k) - sn*s(k)
+ g = sn*s(k+1)
+ s(k+1) = cs*s(k+1)
+ if (wantv) call drot(p,v(1,k),1,v(1,k+1),1,cs,sn)
+ call drotg(f,g,cs,sn)
+ s(k) = f
+ f = cs*e(k) + sn*s(k+1)
+ s(k+1) = -sn*e(k) + cs*s(k+1)
+ g = sn*e(k+1)
+ e(k+1) = cs*e(k+1)
+ if (wantu .and. k .lt. n)
+ * call drot(n,u(1,k),1,u(1,k+1),1,cs,sn)
+ 560 continue
+ e(m-1) = f
+ iter = iter + 1
+ go to 610
+c
+c convergence.
+c
+ 570 continue
+c
+c make the singular value positive.
+c
+ if (s(l) .ge. 0.0d0) go to 580
+ s(l) = -s(l)
+ if (wantv) call dscal(p,-1.0d0,v(1,l),1)
+ 580 continue
+c
+c order the singular value.
+c
+ 590 if (l .eq. mm) go to 600
+c ...exit
+ if (s(l) .ge. s(l+1)) go to 600
+ t = s(l)
+ s(l) = s(l+1)
+ s(l+1) = t
+ if (wantv .and. l .lt. p)
+ * call dswap(p,v(1,l),1,v(1,l+1),1)
+ if (wantu .and. l .lt. n)
+ * call dswap(n,u(1,l),1,u(1,l+1),1)
+ l = l + 1
+ go to 590
+ 600 continue
+ iter = 0
+ m = m - 1
+ 610 continue
+ go to 360
+ 620 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dtrco.f b/com.oracle.truffle.r.native/src/dtrco.f
new file mode 100644
index 0000000000000000000000000000000000000000..8f693e6a818fc23902117385b68d651a726fb556
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dtrco.f
@@ -0,0 +1,161 @@
+ subroutine dtrco(t,ldt,n,rcond,z,job)
+ integer ldt,n,job
+ double precision t(ldt,*),z(*)
+ double precision rcond
+c
+c dtrco estimates the condition of a double precision triangular
+c matrix.
+c
+c on entry
+c
+c t double precision(ldt,n)
+c t contains the triangular matrix. the zero
+c elements of the matrix are not referenced, and
+c the corresponding elements of the array can be
+c used to store other information.
+c
+c ldt integer
+c ldt is the leading dimension of the array t.
+c
+c n integer
+c n is the order of the system.
+c
+c job integer
+c = 0 t is lower triangular.
+c = nonzero t is upper triangular.
+c
+c on return
+c
+c rcond double precision
+c an estimate of the reciprocal condition of t .
+c for the system t*x = b , relative perturbations
+c in t and b of size epsilon may cause
+c relative perturbations in x of size epsilon/rcond .
+c if rcond is so small that the logical expression
+c 1.0 + rcond .eq. 1.0
+c is true, then t may be singular to working
+c precision. in particular, rcond is zero if
+c exact singularity is detected or the estimate
+c underflows.
+c
+c z double precision(n)
+c a work vector whose contents are usually unimportant.
+c if t is close to a singular matrix, then z is
+c an approximate null vector in the sense that
+c norm(a*z) = rcond*norm(a)*norm(z) .
+c
+c linpack. this version dated 08/14/78 .
+c cleve moler, university of new mexico, argonne national lab.
+c
+c subroutines and functions
+c
+c blas daxpy,dscal,dasum
+c fortran dabs,dmax1,dsign
+c
+c internal variables
+c
+ double precision w,wk,wkm,ek
+ double precision tnorm,ynorm,s,sm,dasum
+ integer i1,j,j1,j2,k,kk,l
+ logical lower
+c
+ lower = job .eq. 0
+c
+c compute 1-norm of t
+c
+ tnorm = 0.0d0
+ do 10 j = 1, n
+ l = j
+ if (lower) l = n + 1 - j
+ i1 = 1
+ if (lower) i1 = j
+ tnorm = dmax1(tnorm,dasum(l,t(i1,j),1))
+ 10 continue
+c
+c rcond = 1/(norm(t)*(estimate of norm(inverse(t)))) .
+c estimate = norm(z)/norm(y) where t*z = y and trans(t)*y = e .
+c trans(t) is the transpose of t .
+c the components of e are chosen to cause maximum local
+c growth in the elements of y .
+c the vectors are frequently rescaled to avoid overflow.
+c
+c solve trans(t)*y = e
+c
+ ek = 1.0d0
+ do 20 j = 1, n
+ z(j) = 0.0d0
+ 20 continue
+ do 100 kk = 1, n
+ k = kk
+ if (lower) k = n + 1 - kk
+ if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
+ if (dabs(ek-z(k)) .le. dabs(t(k,k))) go to 30
+ s = dabs(t(k,k))/dabs(ek-z(k))
+ call dscal(n,s,z,1)
+ ek = s*ek
+ 30 continue
+ wk = ek - z(k)
+ wkm = -ek - z(k)
+ s = dabs(wk)
+ sm = dabs(wkm)
+ if (t(k,k) .eq. 0.0d0) go to 40
+ wk = wk/t(k,k)
+ wkm = wkm/t(k,k)
+ go to 50
+ 40 continue
+ wk = 1.0d0
+ wkm = 1.0d0
+ 50 continue
+ if (kk .eq. n) go to 90
+ j1 = k + 1
+ if (lower) j1 = 1
+ j2 = n
+ if (lower) j2 = k - 1
+ do 60 j = j1, j2
+ sm = sm + dabs(z(j)+wkm*t(k,j))
+ z(j) = z(j) + wk*t(k,j)
+ s = s + dabs(z(j))
+ 60 continue
+ if (s .ge. sm) go to 80
+ w = wkm - wk
+ wk = wkm
+ do 70 j = j1, j2
+ z(j) = z(j) + w*t(k,j)
+ 70 continue
+ 80 continue
+ 90 continue
+ z(k) = wk
+ 100 continue
+ s = 1.0d0/dasum(n,z,1)
+ call dscal(n,s,z,1)
+c
+ ynorm = 1.0d0
+c
+c solve t*z = y
+c
+ do 130 kk = 1, n
+ k = n + 1 - kk
+ if (lower) k = kk
+ if (dabs(z(k)) .le. dabs(t(k,k))) go to 110
+ s = dabs(t(k,k))/dabs(z(k))
+ call dscal(n,s,z,1)
+ ynorm = s*ynorm
+ 110 continue
+ if (t(k,k) .ne. 0.0d0) z(k) = z(k)/t(k,k)
+ if (t(k,k) .eq. 0.0d0) z(k) = 1.0d0
+ i1 = 1
+ if (lower) i1 = k + 1
+ if (kk .ge. n) go to 120
+ w = -z(k)
+ call daxpy(n-kk,w,t(i1,k),1,z(i1),1)
+ 120 continue
+ 130 continue
+c make znorm = 1.0
+ s = 1.0d0/dasum(n,z,1)
+ call dscal(n,s,z,1)
+ ynorm = s*ynorm
+c
+ if (tnorm .ne. 0.0d0) rcond = ynorm/tnorm
+ if (tnorm .eq. 0.0d0) rcond = 0.0d0
+ return
+ end
diff --git a/com.oracle.truffle.r.native/src/dtrsl.f b/com.oracle.truffle.r.native/src/dtrsl.f
new file mode 100644
index 0000000000000000000000000000000000000000..b656627dfd5640eceb310ac6d8bffc1a209cdfac
--- /dev/null
+++ b/com.oracle.truffle.r.native/src/dtrsl.f
@@ -0,0 +1,141 @@
+c Triangular Solve dtrsl()
+c ----------------
+c solves systems of the form
+c
+c t * x = b
+c or
+c trans(t) * x = b
+c
+c where t is a triangular matrix of order n. here trans(t)
+c denotes the transpose of the matrix t.
+c
+c on entry
+c
+c t double precision(ldt,n)
+c t contains the matrix of the system. the zero
+c elements of the matrix are not referenced, and
+c the corresponding elements of the array can be
+c used to store other information.
+c
+c ldt integer
+c ldt is the leading dimension of the array t.
+c
+c n integer
+c n is the order of the system.
+c
+c b double precision(n).
+c b contains the right hand side of the system.
+c
+c job integer
+c job specifies what kind of system is to be solved.
+c if job is
+c
+c 00 solve t*x=b, t lower triangular,
+c 01 solve t*x=b, t upper triangular,
+c 10 solve trans(t)*x=b, t lower triangular,
+c 11 solve trans(t)*x=b, t upper triangular.
+c
+c on return
+c
+c b b contains the solution, if info .eq. 0.
+c otherwise b is unaltered.
+c
+c info integer
+c info contains zero if the system is nonsingular.
+c otherwise info contains the index of
+c the first zero diagonal element of t.
+c
+c linpack. this version dated 08/14/78 .
+c g. w. stewart, university of maryland, argonne national lab.
+c
+c subroutines and functions
+c
+c blas: daxpy,ddot
+c fortran mod
+c
+ subroutine dtrsl(t,ldt,n,b,job,info)
+ integer ldt,n,job,info
+ double precision t(ldt,*),b(*)
+c
+c internal variables
+c
+ double precision ddot,temp
+ integer case,j,jj
+c
+c begin block permitting ...exits to 150
+c
+c check for zero diagonal elements.
+c
+ do 10 info = 1, n
+ if (t(info,info) .eq. 0.0d0) go to 150
+c ......exit
+ 10 continue
+ info = 0
+c
+c determine the task and go to it.
+c
+ case = 1
+ if (mod(job,10) .ne. 0) case = 2
+ if (mod(job,100)/10 .ne. 0) case = case + 2
+ go to (20,50,80,110), case
+c
+C Case 1 (job = 00):
+c solve t*x=b for t lower triangular
+c
+ 20 continue
+ b(1) = b(1)/t(1,1)
+ if (n .ge. 2) then
+ do 30 j = 2, n
+ temp = -b(j-1)
+ call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1)
+ b(j) = b(j)/t(j,j)
+ 30 continue
+ endif
+ go to 140
+c
+C Case 2 (job = 01):
+c solve t*x=b for t upper triangular.
+c
+ 50 continue
+ b(n) = b(n)/t(n,n)
+ if (n .ge. 2) then
+ do 60 jj = 2, n
+ j = n - jj + 1
+ temp = -b(j+1)
+ call daxpy(j,temp,t(1,j+1),1,b(1),1)
+ b(j) = b(j)/t(j,j)
+ 60 continue
+ endif
+ go to 140
+c
+C Case 3 (job = 10):
+c solve trans(t)*x=b for t lower triangular.
+c
+ 80 continue
+ b(n) = b(n)/t(n,n)
+ if (n .ge. 2) then
+ do 90 jj = 2, n
+ j = n - jj + 1
+ b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1)
+ b(j) = b(j)/t(j,j)
+ 90 continue
+ endif
+ go to 140
+c
+C Case 4 (job = 11):
+c solve trans(t)*x=b for t upper triangular.
+c
+ 110 continue
+ b(1) = b(1)/t(1,1)
+ if (n .ge. 2) then
+ do 120 j = 2, n
+ b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1)
+ b(j) = b(j)/t(j,j)
+ 120 continue
+ endif
+C
+ 140 continue
+c EXIT:
+ 150 continue
+ return
+ end
diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java
index ebaf7a689fa043a8a7210e0d519a9582b81fa20a..564c7ca42c231db33b64f823b60dfbdce13e495d 100644
--- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java
+++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java
@@ -396,7 +396,9 @@ public class JNR_RFFIFactory extends RFFIFactory implements RFFI, BaseRFFI, RDer
@SlowPath
private static Linpack createAndLoadLib() {
- return LibraryLoader.create(Linpack.class).load("R");
+ // need to load blas lib as Fortran functions in RDerived lib need it
+ LibraryLoader.create(Linpack.class).load("Rblas");
+ return LibraryLoader.create(Linpack.class).load("RDerived");
}
static Linpack linpack() {
diff --git a/mx.fastr/mx_fastr.py b/mx.fastr/mx_fastr.py
index 52ba5c5f96ccf4bf785787915448859f73d62045..655c980a9c1a7b2c7fb0638a7138927f2fd8f959 100644
--- a/mx.fastr/mx_fastr.py
+++ b/mx.fastr/mx_fastr.py
@@ -257,7 +257,7 @@ def _bench_harness_body(args, vmArgs):
'shootout.knucleotide', 'shootout.mandelbrot-ascii', 'shootout.nbody', 'shootout.pidigits',
'shootout.regexdna', 'shootout.reversecomplement', 'shootout.spectralnorm',
'b25.bench.prog-1', 'b25.bench.prog-2', 'b25.bench.prog-3', 'b25.bench.prog-4', 'b25.bench.prog-5',
- 'b25.bench.matcal-1', 'b25.bench.matcal-2', 'b25.bench.matcal-3', 'b25.bench.matcal-4',
+ 'b25.bench.matcal-1', 'b25.bench.matcal-2', 'b25.bench.matcal-3', 'b25.bench.matcal-4', 'b25.bench.matcal-5',
'b25.bench.matfunc-1', 'b25.bench.matfunc-2', 'b25.bench.matfunc-3', 'b25.bench.matfunc-4']
if vmArgs:
marks = ['--J', vmArgs] + marks