view release on metacpan or  search on metacpan

src/cluster.c  view on Meta::CPAN

    free(mask);
    free(data);
}

/* ---------------------------------------------------------------------- */

static double
find_closest_pair(int n, double** distmatrix, int* ip, int* jp)
/*
This function searches the distance matrix to find the pair with the shortest
distance between them. The indices of the pair are returned in ip and jp; the
distance itself is returned by the function.

n          (input) int
The number of elements in the distance matrix.

distmatrix (input) double**
A ragged array containing the distance matrix. The number of columns in each
row is one less than the row index.

ip         (output) int*
A pointer to the integer that is to receive the first index of the pair with
the shortest distance.

jp         (output) int*
A pointer to the integer that is to receive the second index of the pair with
the shortest distance.
*/
{
    int i, j;
    double temp;
    double distance = distmatrix[1][0];

    *ip = 1;
    *jp = 0;
    for (i = 1; i < n; i++) {
        for (j = 0; j < i; j++) {
            temp = distmatrix[i][j];
            if (temp<distance) {
                distance = temp;
                *ip = i;
                *jp = j;
            }
        }
    }
    return distance;
}

/* ********************************************************************* */

static int
svd(int m, int n, double** u, double w[], double** vt)
/*
 * This subroutine is a translation of the Algol procedure svd, described in
 * G.H. Golub and C. Reinsch:
 * "Singular Value Decomposition and Least Squares Solutions",
 * published in
 * Numerische Mathematik 14(5): page 403-420 (1970)
 * and
 * Handbook for Automatic Computation, Volume II: Linear Algebra,
 * (John H. Wilkinson and C. Reinsch; edited by Friedrich L. Bauer and
 * Alston S. Householder): page 134-151 (1971).
 *                                                                  t
 * This subroutine determines the singular value decomposition A=usv  of a
 * real m by n rectangular matrix, where m is greater * than or equal to n.
 * Householder bidiagonalization and a variant of the QR algorithm are used.
 *
 * On input:
 *
 * m  is the number of rows of A (and u).
 *
 * n  is the number of columns of A (and u) and the order of v.
 *
 * u  contains the rectangular input matrix A to be decomposed.
 *
 * On output:
 *
 * the routine returns an integer ierr equal to
 *  0:  to indicate a normal return,
 *  k:  if the k-th singular value has not been determined after 30 iterations,
 * -1:  if memory allocation fails.
 *
 * w  contains the n (non-negative) singular values of a (the
 *    diagonal elements of s). they are unordered.
 *    If an error exit is made, the singular values should be correct for
 *    indices ierr+1, ierr+2, ..., n.
 *
 * u  contains the matrix u (orthogonal column vectors) of the decomposition.
 *    If an error exit is made, the columns of u corresponding to indices of
 *    correct singular values should be correct.
 *                         t
 * vt contains the matrix v (orthogonal) of the decomposition.
 *    If an error exit is made, the columns of v corresponding to indices of
 *    correct singular values should be correct.
 *
 *
 * Questions and comments should be directed to B. S. Garbow,
 * Applied Mathematics division, Argonne National Laboratory
 *
 * Modified to eliminate machep
 *
 * Translated to C by Michiel de Hoon, Human Genome Center,
 * University of Tokyo, for inclusion in the C Clustering Library.
 * This routine is less general than the original svd routine, as
 * it focuses on the singular value decomposition as needed for
 * clustering. In particular,
 *  - We calculate both u and v in all cases
 *  - We pass the input array A via u; this array is subsequently overwritten.
 *  - We allocate for the array rv1, used as a working space, internally in
 *    this routine, instead of passing it as an argument.
 *    If the allocation fails, svd returns -1.
 * 2003.06.05
 */
{
    int i, j, k, i1, k1, l1, its;
    double c, f, h, s, x, y, z;
    int l = 0;
    int ierr = 0;
    double g = 0.0;
    double scale = 0.0;
    double anorm = 0.0;