package com.github.fommil.netlib;
import com.google.common.base.Stopwatch;
import lombok.extern.java.Log;
import java.util.concurrent.TimeUnit;
/**
* Routines that give a bit of a workout of BLAS.
*
* @author Sam Halliday (adapted to netlib-java)
* @author David M. Doolin (bugfixes)
* @author Jonathan Hardwick (Java optimisations)
* @author Reed Wade (Java translation)
* @author Jack Dongarra (bugfixes)
* @author Bonnie Toy (C translation)
* @see <a href="http://www.netlib.org/linpack/">LINPACK</a>
*/
@Log
public class Linpack implements Benchmark {
private BLAS blas = BLAS.getInstance();
public static void main(String[] args) {
Linpack linpack = new Linpack();
linpack.benchmark();
}
@Override
public long benchmark() {
double a[][] = new double[200][201];
double b[] = new double[200];
int n = 100, lda = 201;
int ipvt[] = new int[200];
double ops = (2.0 * (n * n * n)) / 3.0 + 2.0 * (n * n);
matgen(a, n, b);
Stopwatch watch = new Stopwatch();
watch.start();
dgefa(a, lda, n, ipvt);
dgesl(a, lda, n, ipvt, b, 0);
watch.stop();
long total = watch.elapsed(TimeUnit.NANOSECONDS);
double mflops = 1000 * ops / total;
log.info("Mflops: " + mflops);
return total;
}
final double matgen(double a[][], int n, double b[]) {
double norma;
int init, i, j;
init = 1325;
norma = 0.0;
/* Next two for() statements switched. Solver wants
matrix in column order. --dmd 3/3/97
*/
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
init = 3125 * init % 65536;
a[j][i] = (init - 32768.0) / 16384.0;
norma = (a[j][i] > norma) ? a[j][i] : norma;
}
}
for (i = 0; i < n; i++) {
b[i] = 0.0;
}
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
b[i] += a[j][i];
}
}
return norma;
}
/*
dgefa factors a double precision matrix by gaussian elimination.
dgefa is usually called by dgeco, but it can be called
directly with a saving in time if rcond is not needed.
(time for dgeco) = (1 + 9/n)*(time for dgefa) .
on entry
a double precision[n][lda]
the matrix to be factored.
lda integer
the leading dimension of the array a .
n integer
the order of the matrix a .
on return
a an upper triangular matrix and the multipliers
which were used to obtain it.
the factorization can be written a = l*u where
l is a product of permutation and unit lower
triangular matrices and u is upper triangular.
ipvt integer[n]
an integer vector of pivot indices.
info integer
= 0 normal value.
= k if u[k][k] .eq. 0.0 . this is not an error
condition for this subroutine, but it does
indicate that dgesl or dgedi will divide by zero
if called. use rcond in dgeco for a reliable
indication of singularity.
linpack. this version dated 08/14/78.
cleve moler, university of new mexico, argonne national lab.
functions
blas daxpy,dscal,idamax
*/
final int dgefa(double a[][], int lda, int n, int ipvt[]) {
double[] col_k, col_j;
double t;
int j, k, kp1, l, nm1;
int info;
// gaussian elimination with partial pivoting
info = 0;
nm1 = n - 1;
if (nm1 >= 0) {
for (k = 0; k < nm1; k++) {
col_k = a[k];
kp1 = k + 1;
// find l = pivot index
l = blas.idamax(n - k, col_k, k, 1) + k;
ipvt[k] = l;
// zero pivot implies this column already triangularized
if (col_k[l] != 0) {
// interchange if necessary
if (l != k) {
t = col_k[l];
col_k[l] = col_k[k];
col_k[k] = t;
}
// compute multipliers
t = -1.0 / col_k[k];
blas.dscal(n - (kp1), t, col_k, kp1, 1);
// row elimination with column indexing
for (j = kp1; j < n; j++) {
col_j = a[j];
t = col_j[l];
if (l != k) {
col_j[l] = col_j[k];
col_j[k] = t;
}
blas.daxpy(n - (kp1), t, col_k, kp1, 1,
col_j, kp1, 1);
}
} else {
info = k;
}
}
}
ipvt[n - 1] = n - 1;
if (a[(n - 1)][(n - 1)] == 0) info = n - 1;
return info;
}
/*
dgesl solves the double precision system
a * x = b or trans(a) * x = b
using the factors computed by dgeco or dgefa.
on entry
a double precision[n][lda]
the output from dgeco or dgefa.
lda integer
the leading dimension of the array a .
n integer
the order of the matrix a .
ipvt integer[n]
the pivot vector from dgeco or dgefa.
b double precision[n]
the right hand side vector.
job integer
= 0 to solve a*x = b ,
= nonzero to solve trans(a)*x = b where
trans(a) is the transpose.
on return
b the solution vector x .
error condition
a division by zero will occur if the input factor contains a
zero on the diagonal. technically this indicates singularity
but it is often caused by improper arguments or improper
setting of lda . it will not occur if the subroutines are
called correctly and if dgeco has set rcond .gt. 0.0
or dgefa has set info .eq. 0 .
to compute inverse(a) * c where c is a matrix
with p columns
dgeco(a,lda,n,ipvt,rcond,z)
if (!rcond is too small){
for (j=0,j<p,j++)
dgesl(a,lda,n,ipvt,c[j][0],0);
}
linpack. this version dated 08/14/78 .
cleve moler, university of new mexico, argonne national lab.
functions
blas daxpy,ddot
*/
final void dgesl(double a[][], int lda, int n, int ipvt[], double b[], int job) {
double t;
int k, kb, l, nm1, kp1;
nm1 = n - 1;
if (job == 0) {
// job = 0 , solve a * x = b. first solve l*y = b
if (nm1 >= 1) {
for (k = 0; k < nm1; k++) {
l = ipvt[k];
t = b[l];
if (l != k) {
b[l] = b[k];
b[k] = t;
}
kp1 = k + 1;
blas.daxpy(n - (kp1), t, a[k], kp1, 1, b, kp1, 1);
}
}
// now solve u*x = y
for (kb = 0; kb < n; kb++) {
k = n - (kb + 1);
b[k] /= a[k][k];
t = -b[k];
blas.daxpy(k, t, a[k], 0, 1, b, 0, 1);
}
} else {
// job = nonzero, solve trans(a) * x = b. first solve trans(u)*y = b
for (k = 0; k < n; k++) {
t = blas.ddot(k, a[k], 0, 1, b, 0, 1);
b[k] = (b[k] - t) / a[k][k];
}
// now solve trans(l)*x = y
if (nm1 >= 1) {
for (kb = 1; kb < nm1; kb++) {
k = n - (kb + 1);
kp1 = k + 1;
b[k] += blas.ddot(n - (kp1), a[k], kp1, 1, b, kp1, 1);
l = ipvt[k];
if (l != k) {
t = b[l];
b[l] = b[k];
b[k] = t;
}
}
}
}
}
}