Dlasrt.java
package org.mklab.sdpj.gpack.lapackwrap;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.gpack.blaswrap.BLAS;
/**
* @author koga
* @version $Revision$, 2009/04/25
* @param <RS> type of real scalar
* @param <RM> type of real matrix
* @param <CS> type of complex scalar
* @param <CM> type of complex Matrix
*/
public class Dlasrt<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {
/* -- LAPACK routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
September 30, 1994
Purpose
=======
Sort the numbers in D in increasing order (if ID = 'I') or
in decreasing order (if ID = 'D' ).
Use Quick Sort, reverting to Insertion sort on arrays of
size <= 20. Dimension of STACK limits N to about 2**32.
Arguments
=========
ID (input) CHARACTER*1
= 'I': sort D in increasing order;
= 'D': sort D in decreasing order.
N (input) INTEGER
The length of the array D.
D (input/output) DOUBLE PRECISION array, dimension (N)
On entry, the array to be sorted.
On exit, D has been sorted into increasing order
(D(1) <= ... <= D(N) ) or into decreasing order
(D(1) >= ... >= D(N) ), depending on ID.
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
=====================================================================
*/
/**
* @param id input
* @param n input
* @param d input/output
* @return result
*/
int execute(String id, int n, RS[] d) {
int[] stack = new int[64]; /* was [2][32] */
int pointer_d = -1;
int pointer_stack = 0;
int info = 0;
int dir = -1;
if (BLAS.lsame(id, "D")) { //$NON-NLS-1$
dir = 0;
} else if (BLAS.lsame(id, "I")) { //$NON-NLS-1$
dir = 1;
}
if (dir == -1) {
info = -1;
} else if (n < 0) {
info = -2;
}
if (info != 0) {
BLAS.xerbla("DLASRT", -(info)); //$NON-NLS-1$
return info;
}
/* Quick return if possible */
if (n <= 1) {
return info;
}
int stkpnt = 1;
stack[2 + 1 - 3 + pointer_stack] = 1;
stack[2 + 2 - 3 + pointer_stack] = n;
do {
int start = stack[(stkpnt) * 2 + 1 - 3 + pointer_stack];
int endd = stack[(stkpnt) * 2 + 2 - 3 + pointer_stack];
--stkpnt;
if (endd - start <= 20 && endd - start > 0) {
/* Do Insertion sort on D( START:ENDD ) */
if (dir == 0) {
/* Sort into decreasing order */
for (int i = start + 1; i <= endd; ++i) {
for (int j = i; j >= start + 1; --j) {
if (d[pointer_d + j].isGreaterThan(d[pointer_d + j - 1])) {
RS dmnmx = d[pointer_d + j];
d[pointer_d + j] = d[pointer_d + j - 1];
d[pointer_d + j - 1] = dmnmx;
} else {
break;
}
/* L20: */
}
}
} else {
/* Sort into increasing order */
for (int i = start + 1; i <= endd; ++i) {
for (int j = i; j >= start + 1; --j) {
if (d[pointer_d + j].isLessThan(d[pointer_d + j - 1])) {
RS dmnmx = d[pointer_d + j];
d[pointer_d + j] = d[pointer_d + j - 1];
d[pointer_d + j - 1] = dmnmx;
} else {
break;
}
/* L40: */
}
}
}
} else if (endd - start > 20) {
/* Partition D( START:ENDD ) and stack parts, largest one first
Choose partition entry as median of 3 */
RS d1 = d[pointer_d + start];
RS d2 = d[pointer_d + endd];
RS d3 = d[pointer_d + (start + endd) / 2];
RS dmnmx;
if (d1.isLessThan(d2)) {
if (d3.isLessThan(d1)) {
dmnmx = d1;
} else if (d3.isLessThan(d2)) {
dmnmx = d3;
} else {
dmnmx = d2;
}
} else {
if (d3.isLessThan(d2)) {
dmnmx = d2;
} else if (d3.isLessThan(d1)) {
dmnmx = d3;
} else {
dmnmx = d1;
}
}
if (dir == 0) {
/* Sort into decreasing order */
int i = start - 1;
int j = endd + 1;
//L60:
do {
//L70:
do {
--j;
} while (d[pointer_d + j].isLessThan(dmnmx));
//L80:
do {
++i;
} while (d[pointer_d + i].isGreaterThan(dmnmx));
if (i < j) {
RS tmp = d[pointer_d + i];
d[pointer_d + i] = d[pointer_d + j];
d[pointer_d + j] = tmp;
}
} while (i < j);
if (j - start > endd - j - 1) {
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = start;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = j;
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = j + 1;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = endd;
} else {
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = j + 1;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = endd;
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = start;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = j;
}
} else {
/* Sort into increasing order */
int i = start - 1;
int j = endd + 1;
//L90:
do {
//L100:
do {
--j;
} while (d[pointer_d + j].isGreaterThan(dmnmx));
//L110:
do {
++i;
} while (d[pointer_d + i].isLessThan(dmnmx));
if (i < j) {
RS tmp = d[pointer_d + i];
d[pointer_d + i] = d[pointer_d + j];
d[pointer_d + j] = tmp;
}
} while (i < j);
if (j - start > endd - j - 1) {
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = start;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = j;
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = j + 1;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = endd;
} else {
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = j + 1;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = endd;
++stkpnt;
stack[(stkpnt) * 2 + 1 - 3 + pointer_stack] = start;
stack[(stkpnt) * 2 + 2 - 3 + pointer_stack] = j;
}
}
}
} while (stkpnt > 0);
return info;
}
}