NAME

ICC::Support::Lapack - Perl interface to LAPACK and BLAS libraries, used by ICC::Profile modules.

SYNOPSIS

    use ICC::Support::Lapack;

    # interface to selected LAPACK routines
    $info = ICC::Support::Lapack::dgels($trans, $m, $n, $nrhs, $a, $lda, $b, $ldb);
    $info = ICC::Support::Lapack::dgesv($n, $nrhs, $a, $lda, $b, $ldb);
    $info = ICC::Support::Lapack::dgesvd($jobu, $jobvt, $m, $n, $a, $lda, $s, $u, $ldu, $vt, $ldvt);
    $info = ICC::Support::Lapack::dgeqrf($m, $n, $a, $lda, $tau);

    # interface to selected BLAS routines
    ICC::Support::Lapack::dgemm($transa, $transb, $m, $n, $k, $alpha, $a, $lda, $b, $ldb, $beta, $c, $ldc);
    ICC::Support::Lapack::dgemv($trans, $m, $n, $alpha, $a, $lda, $x, $incx, $beta, $y, $incy);

    # functions used by clut.pm module
    $y = ICC::Support::Lapack::clut_mat_trans($gsa, $x, $cache);

    # functions used by matf.pm module
    ($info, $ab) = ICC::Support::Lapack::matf_fit($x, $y, $offset);
    $y = ICC::Support::Lapack::matf_mat_trans($x, $a, [$b]);
    $y = ICC::Support::Lapack::matf_vec_trans($x, $a, [$b]);
    ($info, $x) = ICC::Support::Lapack::matf_inv($y, $a, [$b]);

    # functions used by nMIX.pm module
    $y = ICC::Support::Lapack::nMIX_mat_trans($x, $cache, [$delta]);
    $y = ICC::Support::Lapack::nMIX_vec_trans($x, $cache, [$delta]);
    $y = ICC::Support::Lapack::nMIX_jacobian($x, $cache, [$delta, $out]);
    $y = ICC::Support::Lapack::nMIX_power($x, $delta);
    $cache = ICC::Support::Lapack::cache_2D($coef);

    # functions used by nNET.pm module
    ($info, $ab) = ICC::Support::Lapack::nNET_fit($x, $y, $offset);

    # functions used by nPINT.pm module
    ($info, $coef) = ICC::Support::Lapack::nPINT_fit($pda, $x, $y);
    $y = ICC::Support::Lapack::nPINT_mat_trans($pda, $x, $cache);
    $y = ICC::Support::Lapack::nPINT_vec_trans($pda, $x, $cache);
    $jac = ICC::Support::Lapack::nPINT_jacobian($pda, $x, $cache);
    $cols = ICC::Support::Lapack::xcols($pda);
    $cache = ICC::Support::Lapack::cache_2D($coef);

    # matrix functions
    $zeros = ICC::Support::Lapack::zeros($m, [$n]);
    $ones = ICC::Support::Lapack::ones($m, [$n]);
    $idn = ICC::Support::Lapack::eye($m, [$n]);
    $y = ICC::Support::Lapack::power($a, $v);
    ($info, $x) = ICC::Support::Lapack::chol($a);
    ($info, $x) = ICC::Support::Lapack::solve($a, $b);
    ($info, $x) = ICC::Support::Lapack::trisolve($dl, $d, $du, $b);
    ($info, $x) = ICC::Support::Lapack::inv($a);
    ($info, $delta) = ICC::Support::Lapack::normal($a, $b, [$w]);
    ($info, $q, $r) = ICC::Support::Lapack::qr($a, [$eco]);
    $y = ICC::Support::Lapack::dot($a, $b);
    $y = ICC::Support::Lapack::vec_xplus($a, $x, [$c]);
    $y = ICC::Support::Lapack::mat_xplus($a, $b, [$c]);

    # statistics functions
    $mean = ICC::Support::Lapack::mean($a);
    ($cov, $mean) = ICC::Support::Lapack::cov($a);
    ($info, $u, $sigma, $vt, $mean) = ICC::Support::Lapack::pca($a);
    $dist = ICC::Support::Lapack::mahal($v1, $v2, $chol);
    $dist = ICC::Support::Lapack::euclid($v1, $v2);
    $sn = ICC::Support::Lapack::skew_normal($x, $location, $scale, $shape);

    # special functions
    $exp = ICC::Support::Lapack::expm1($x);
    $log = ICC::Support::Lapack::log1p($x);
    $erf = ICC::Support::Lapack::erf($x);
    $gamma = ICC::Support::Lapack::tgamma($x);
    $beta = ICC::Support::Lapack::beta($x, $y);

    # miscellaneous functions
    $round = ICC::Support::Lapack::round($x);
    $y = ICC::Support::Lapack::bandpass($x, $mu, $sig1, [$sig2, [$gam1, [$gam2]]]);
    $y = ICC::Support::Lapack::hermite($t, $y0, $y1, $m0, $m1);
    $y = ICC::Support::Lapack::bernstein($t, $coef);

DESCRIPTION

ICC::Support::Lapack is a Perl interface to the LAPACK and BLAS libraries, used by ICC::Profile modules.

LAPACK (Linear Algebra PACKage) provides routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems. The associated matrix factorizations (LU, Cholesky, QR, SVD, Schur, generalized Schur) are also provided.

The BLAS (Basic Linear Algebra Subprograms) are routines that provide standard building blocks for performing basic vector and matrix operations. The Level 1 BLAS perform scalar, vector and vector-vector operations, the Level 2 BLAS perform matrix-vector operations, and the Level 3 BLAS perform matrix-matrix operations.

LAPACK and BLAS are written in FORTRAN and provide routines for real and complex matrices, in single and double precision. Currently, LAPACK has 1638 routines and BLAS has 151 routines. This module only interfaces to the routines used by the ICC::Profile modules. If you need an interface to other routines, it would be easy to modify Lapack.xs using the existing code as a model.

This module also includes functions specific to the clut.pm, matf.pm, nMIX.pm, nNET.pm, and nPINT.pm modules, and some matrix and statistical functions, similar to those in MATLAB or R.

FUNCTIONS

Interface to Selected LAPACK Routines

dgels

DGELS solves overdetermined or underdetermined real linear systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A. It is assumed that A has full rank.

The following options are provided:

    1. If TRANS = 'N' and m >= n:  find the least squares solution of
       an overdetermined system, i.e., solve the least squares problem
                    minimize || B - A*X ||.

    2. If TRANS = 'N' and m < n:  find the minimum norm solution of
       an underdetermined system A * X = B.

    3. If TRANS = 'T' and m >= n:  find the minimum norm solution of
       an undetermined system A**T * X = B.

    4. If TRANS = 'T' and m < n:  find the least squares solution of
       an overdetermined system, i.e., solve the least squares problem
                    minimize || B - A**T * X ||.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix X.

Arguments

    TRANS   (input) CHARACTER*1
            = 'N': the linear system involves A;
            = 'T': the linear system involves A**T.

    M       (input) INTEGER
            The number of rows of the matrix A. M >= 0.

    N       (input) INTEGER
            The number of columns of the matrix A. N >= 0.

    NRHS    (input) INTEGER
            The number of right hand sides, i.e., the number of
            columns of the matrices B and X. NRHS >=0.

    A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
            On entry, the M-by-N matrix A.
            On exit,
              if M >= N, A is overwritten by details of its QR
                         factorization as returned by DGEQRF;
              if M <  N, A is overwritten by details of its LQ
                         factorization as returned by DGELQF.

    LDA     (input) INTEGER
            The leading dimension of the array A. LDA >= max(1,M).

    B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
            On entry, the matrix B of right hand side vectors, stored
            columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
            if TRANS = 'T'.
            On exit, if INFO = 0, B is overwritten by the solution
            vectors, stored columnwise:
            if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
            squares solution vectors; the residual sum of squares for the
            solution in each column is given by the sum of squares of
            elements N+1 to M in that column;
            if TRANS = 'N' and m < n, rows 1 to N of B contain the
            minimum norm solution vectors;
            if TRANS = 'T' and m >= n, rows 1 to M of B contain the
            minimum norm solution vectors;
            if TRANS = 'T' and m < n, rows 1 to M of B contain the
            least squares solution vectors; the residual sum of squares
            for the solution in each column is given by the sum of
            squares of elements M+1 to N in that column.

    LDB     (input) INTEGER
            The leading dimension of the array B. LDB >= MAX(1,M,N).

    INFO    (output) INTEGER
            = 0:  successful exit
            < 0:  if INFO = -i, the i-th argument had an illegal value
            > 0:  if INFO =  i, the i-th diagonal element of the
                  triangular factor of A is zero, so that A does not have
                  full rank; the least squares solution could not be
                  computed.

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $trans = 'N';
    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1], [-1, 11, 12]];
    $b = [[0], [5], [2], [7]];
    $m = @{$a};
    $n = @{$a->[0]};
    $nrhs = @{$b->[0]};
    $lda = $m;
    $ldb = $m;
    $info = ICC::Support::Lapack::dgels($trans, $m, $n, $nrhs, $a, $lda, $b, $ldb);
    print Dumper($info, $b);

dgesv

DGESV computes the solution to a real system of linear equations

    A * X = B

where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

The LU decomposition with partial pivoting and row interchanges is used to factor A as

    A = P * L * U

where P is a permutation matrix, L is unit lower triangular, and U is upper triangular. The factored form of A is then used to solve the system of equations A * X = B.

Arguments

    N       (input) INTEGER
            The number of linear equations, i.e., the order of the
            matrix A. N >= 0.

    NRHS    (input) INTEGER
            The number of right hand sides, i.e., the number of columns
            of the matrix B. NRHS >= 0.

    A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
            On entry, the N-by-N coefficient matrix A.
            On exit, the factors L and U from the factorization
            A = P*L*U; the unit diagonal elements of L are not stored.

    LDA     (input) INTEGER
            The leading dimension of the array A. LDA >= max(1,N).

    B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
            On entry, the N-by-NRHS matrix of right hand side matrix B.
            On exit, if INFO = 0, the N-by-NRHS solution matrix X.

    LDB     (input) INTEGER
            The leading dimension of the array B. LDB >= max(1,N).

    INFO    (output) INTEGER
            = 0:  successful exit
            < 0:  if INFO = -i, the i-th argument had an illegal value
            > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
                  has been completed, but the factor U is exactly
                  singular, so the solution could not be computed.

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1]];
    $b = [[0], [5], [2]];
    $n = @{$a};
    $nrhs = @{$b->[0]};
    $lda = $n;
    $ldb = $n;
    $info = ICC::Support::Lapack::dgesv($n, $nrhs, $a, $lda, $b, $ldb);
    print Dumper($info, $b);

dgesvd

DGESVD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written

    A = U * SIGMA * transpose(V)

where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A.

Note that the routine returns V**T, not V.

Arguments

    JOBU    (input) CHARACTER*1
            Specifies options for computing all or part of the matrix U:
            = 'A':  all M columns of U are returned in array U:
            = 'S':  the first min(m,n) columns of U (the left singular
                    vectors) are returned in the array U;
            = 'O':  the first min(m,n) columns of U (the left singular
                    vectors) are overwritten on the array A;
            = 'N':  no columns of U (no left singular vectors) are
                    computed.

    JOBVT   (input) CHARACTER*1
            Specifies options for computing all or part of the matrix
            V**T:
            = 'A':  all N rows of V**T are returned in the array VT;
            = 'S':  the first min(m,n) rows of V**T (the right singular
                    vectors) are returned in the array VT;
            = 'O':  the first min(m,n) rows of V**T (the right singular
                    vectors) are overwritten on the array A;
            = 'N':  no rows of V**T (no right singular vectors) are
                    computed.

            JOBVT and JOBU cannot both be 'O'.

    M       (input) INTEGER
            The number of rows of the input matrix A. M >= 0.

    N       (input) INTEGER
            The number of columns of the input matrix A. N >= 0.

    A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
            On entry, the M-by-N matrix A.
            On exit,
            if JOBU = 'O',  A is overwritten with the first min(m,n)
                            columns of U (the left singular vectors,
                            stored columnwise);
            if JOBVT = 'O', A is overwritten with the first min(m,n)
                            rows of V**T (the right singular vectors,
                            stored rowwise);
            if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
                            are destroyed.

    LDA     (input) INTEGER
            The leading dimension of the array A. LDA >= max(1,M).

    S       (output) DOUBLE PRECISION array, dimension (min(M,N))
            The singular values of A, sorted so that S(i) >= S(i+1).

    U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
            (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
            If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
            if JOBU = 'S', U contains the first min(m,n) columns of U
            (the left singular vectors, stored columnwise);
            if JOBU = 'N' or 'O', U is not referenced.

    LDU     (input) INTEGER
            The leading dimension of the array U. LDU >= 1; if
            JOBU = 'S' or 'A', LDU >= M.

    VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
            If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
            V**T;
            if JOBVT = 'S', VT contains the first min(m,n) rows of
            V**T (the right singular vectors, stored rowwise);
            if JOBVT = 'N' or 'O', VT is not referenced.

    LDVT    (input) INTEGER
            The leading dimension of the array VT. LDVT >= 1; if
            JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).

    INFO    (output) INTEGER
            = 0:  successful exit.
            < 0:  if INFO = -i, the i-th argument had an illegal value.
            > 0:  if DBDSQR did not converge, INFO specifies how many
                  superdiagonals of an intermediate bidiagonal form B
                  did not converge to zero. See the description of WORK
                  above for details.

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $jobu = 'A';
    $jobvt = 'A';
    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1]];
    $m = @{$a};
    $n = @{$a->[0]};
    $lda = $m;
    $s = [];
    $u = [];
    $ldu = $m;
    $vt = [];
    $ldvt = $n;
    $info = ICC::Support::Lapack::dgesvd($jobu, $jobvt, $m, $n, $a, $lda, $s, $u, $ldu, $vt, $ldvt);
    print Dumper($info, $s, $u, $vt);

dgeqrf

DGEQRF computes a QR factorization of a real M-by-N matrix A:

    A = Q * R.

Arguments

    M       (input) INTEGER
            The number of rows of the matrix A. M >= 0.

    N       (input) INTEGER
            The number of columns of the matrix A. N >= 0.

    A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
            On entry, the M-by-N matrix A.
            On exit, the elements on and above the diagonal of the array
            contain the min(M,N)-by-N upper trapezoidal matrix R (R is
            upper triangular if m >= n); the elements below the diagonal,
            with the array TAU, represent the orthogonal matrix Q as a
            product of min(m,n) elementary reflectors (see Further
            Details).

    LDA     (input) INTEGER
            The leading dimension of the array A. LDA >= max(1,M).

    TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
            The scalar factors of the elementary reflectors (see Further
            Details).

    INFO    (output) INTEGER
            = 0:  successful exit
            < 0:  if INFO = -i, the i-th argument had an illegal value

Further Details

The matrix Q is represented as a product of elementary reflectors

    Q = H(1) H(2) . . . H(k), where k = min(m,n).

Each H(i) has the form

    H(i) = I - tau * v * v**T

where tau is a real scalar, and v is a real vector with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), and tau in TAU(i).

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1], [-1, 11, 12]];
    $m = @{$a};
    $n = @{$a->[0]};
    $lda = $m;
    $tau = [];
    $info = ICC::Support::Lapack::dgeqrf($m, $n, $a, $lda, $tau);
    print Dumper($info, $a, $tau);

Interface to Selected BLAS Routines

dgemm

DGEMM performs one of the matrix-matrix operations

    C : = alpha*op( A )*op( B ) + beta*C,

where op( X ) is one of

    op( X ) = X   or   op( X ) = X**T,

alpha and beta are scalars, and A, B and C are matrices, with op( A ) an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.

Arguments

    TRANSA - CHARACTER*1.
             On entry, TRANSA specifies the form of op( A ) to be used in
             the matrix multiplication as follows:

                TRANSA = 'N' or 'n',  op( A ) = A.

                TRANSA = 'T' or 't',  op( A ) = A**T.

                TRANSA = 'C' or 'c',  op( A ) = A**T.

             Unchanged on exit.

    TRANSB - CHARACTER*1.
             On entry, TRANSB specifies the form of op( B ) to be used in
             the matrix multiplication as follows:

                TRANSB = 'N' or 'n',  op( B ) = B.

                TRANSB = 'T' or 't',  op( B ) = B**T.

                TRANSB = 'C' or 'c',  op( B ) = B**T.

             Unchanged on exit.

    M      - INTEGER.
             On entry,  M  specifies  the number  of rows  of the  matrix
             op( A )  and of the  matrix  C. M  must  be at least  zero.
             Unchanged on exit.

    N      - INTEGER.
             On entry,  N  specifies the number  of columns of the matrix
             op( B ) and the number of columns of the matrix C. N must be
             at least zero.
             Unchanged on exit.

    K      - INTEGER.
             On entry,  K  specifies  the number of columns of the matrix
             op( A ) and the number of rows of the matrix op( B ). K must
             be at least  zero.
             Unchanged on exit.

    ALPHA  - DOUBLE PRECISION.
             On entry, ALPHA specifies the scalar alpha.
             Unchanged on exit.

    A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
             k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
             Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
             part of the array  A  must contain the matrix  A,  otherwise
             the leading  k by m  part of the array  A  must contain  the
             matrix A.
             Unchanged on exit.

    LDA    - INTEGER.
             On entry, LDA specifies the first dimension of A as declared
             in the calling (sub) program. When  TRANSA = 'N' or 'n' then
             LDA must be at least  max( 1, m ), otherwise  LDA must be at
             least  max( 1, k ).
             Unchanged on exit.

    B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
             n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
             Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
             part of the array  B  must contain the matrix  B,  otherwise
             the leading  n by k  part of the array  B  must contain  the
             matrix B.
             Unchanged on exit.

    LDB    - INTEGER.
             On entry, LDB specifies the first dimension of B as declared
             in the calling (sub) program. When  TRANSB = 'N' or 'n' then
             LDB must be at least  max( 1, k ), otherwise  LDB must be at
             least  max( 1, n ).
             Unchanged on exit.

    BETA   - DOUBLE PRECISION.
             On entry,  BETA  specifies the scalar  beta. When  BETA  is
             supplied as zero then C need not be set on input.
             Unchanged on exit.

    C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
             Before entry, the leading  m by n  part of the array  C must
             contain the matrix  C,  except when  beta  is zero, in which
             case C need not be set on entry.
             On exit, the array  C  is overwritten by the  m by n  matrix
             ( alpha*op( A )*op( B ) + beta*C ).

    LDC    - INTEGER.
             On entry, LDC specifies the first dimension of C as declared
             in  the  calling  (sub)  program. LDC  must  be  at  least
             max( 1, m ).
             Unchanged on exit.

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $transa = 'N';
    $transb = 'N';
    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1], [-1, 11, 12]];
    $b = [[1, 2, 3, 4], [5, -1, 7, 2], [-1, 8, 9, 3]];
    $c = [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]];
    $m = @{$a};
    $n = @{$b->[0]};
    $k = @{$b};
    $lda = $m;
    $ldb = $k;
    $ldc = $m;
    $alpha = 1;
    $beta = 1;
    ICC::Support::Lapack::dgemm($transa, $transb, $m, $n, $k, $alpha, $a, $lda, $b, $ldb, $beta, $c, $ldc);
    print Dumper($c);

dgemv

DGEMV performs one of the matrix-vector operations

    y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,

where alpha and beta are scalars, x and y are vectors and A is an m by n matrix.

Arguments

    TRANS  - CHARACTER*1.
             On entry, TRANS specifies the operation to be performed as
             follows:

                TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.

                TRANS = 'T' or 't'   y := alpha*A**T*x + beta*y.

                TRANS = 'C' or 'c'   y := alpha*A**T*x + beta*y.

             Unchanged on exit.

    M      - INTEGER.
             On entry, M specifies the number of rows of the matrix A.
             M must be at least zero.
             Unchanged on exit.

    N      - INTEGER.
             On entry, N specifies the number of columns of the matrix A.
             N must be at least zero.
             Unchanged on exit.

    ALPHA  - DOUBLE PRECISION.
             On entry, ALPHA specifies the scalar alpha.
             Unchanged on exit.

    A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
             Before entry, the leading m by n part of the array A must
             contain the matrix of coefficients.
             Unchanged on exit.

    LDA    - INTEGER.
             On entry, LDA specifies the first dimension of A as declared
             in the calling (sub) program. LDA must be at least
             max( 1, m ).
             Unchanged on exit.

    X      - DOUBLE PRECISION array of DIMENSION at least
             ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
             and at least
             ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
             Before entry, the incremented array X must contain the
             vector x.
             Unchanged on exit.

    INCX   - INTEGER.
             On entry, INCX specifies the increment for the elements of
             X. INCX must not be zero.
             Unchanged on exit.

    BETA   - DOUBLE PRECISION.
             On entry, BETA specifies the scalar beta. When BETA is
             supplied as zero then Y need not be set on input.
             Unchanged on exit.

    Y      - DOUBLE PRECISION array of DIMENSION at least
             ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
             and at least
             ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
             Before entry with BETA non-zero, the incremented array Y
             must contain the vector y. On exit, Y is overwritten by the
             updated vector y.

    INCY   - INTEGER.
             On entry, INCY specifies the increment for the elements of
             Y. INCY must not be zero.
             Unchanged on exit.

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $trans = 'N';
    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1], [-1, 11, 12]];
    $x = [1, 2, 3];
    $y = [0, 0, 0, 0];
    $m = @{$a};
    $n = @{$a->[0]};
    $lda = $m;
    $incx = 1;
    $incy = 1;
    $alpha = 1;
    $beta = 1;
    ICC::Support::Lapack::dgemv($trans, $m, $n, $alpha, $a, $lda, $x, $incx, $beta, $y, $incy);
    print Dumper($y);

Functions Used By clut.pm Module

The 'clut.pm' module implements the ICC Profile 'clut' tag. The acronym 'CLUT' stands for 'color lookup table'. A 'clut' object contains a lookup table of variable dimensions and grid size. The 'transform' method locates the grid points surrounding the input value(s), and interpolates the output value(s). The following Lapack functions support the 'transform' and 'jacobian' methods of the 'clut.pm' module.

clut_vec_trans

clut_vec_trans transforms a vector input to a vector output. Uses the BLAS DGEMV function.

    GSA - a N-length vector (grid size array)

    X - a N-length input vector

    CACHE - a P-by-M matrix (CLUT) in cached format (where P is total number of grid points)

    Y - a M-length output vector

Example

    use ICC::Profile;
    use ICC::Support::Lapack;
    use Data::Dumper;

    $profile = ICC::Profile->new('../../test_data/profiles/originals/GRACoL2006_FullK_i1P.icc'); # open a profile
    $A2B0 = $profile->tag('A2B0');
    $clut = $A2B0->clut(); # get clut tag
    $cache = ICC::Support::Lapack::cache_2D($clut->array()); # get cached CLUT
    $gsa = $clut->gsa(); # get gsa

    $x = [0.1, 0.3, 0.5, 0.1]; # input vector
    $y = ICC::Support::Lapack::clut_vec_trans($gsa, $x, $cache); # transform vector
    print Dumper($x, $y);

clut_jacobian

clut_jacobian computes the Jacobian matrix for a single sample. The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMM function.

    GSA - a N-length vector (grid size array)

    X - a N-length input vector

    CACHE - a P-by-M matrix (CLUT) in cached format (where P is total number of grid points)

    J - an M x N Jacobian matrix

Example

    use ICC::Profile;
    use ICC::Support::Lapack;
    use Data::Dumper;

    $profile = ICC::Profile->new('../../test_data/profiles/originals/GRACoL2006_FullK_i1P.icc'); # open a profile
    $A2B0 = $profile->tag('A2B0');
    $clut = $A2B0->clut(); # get clut tag
    $cache = ICC::Support::Lapack::cache_2D($clut->array()); # get cached CLUT
    $gsa = $clut->gsa(); # get gsa

    $x = [0.1, 0.3, 0.5, 0.1]; # input vector
    $jac = ICC::Support::Lapack::clut_jacobian($gsa, $x, $cache); # Jacobian matrix
    print Dumper($x, $jac);

clut_jacobian_ext

clut_jacobian_ext computes the Jacobian matrix for a single sample. Same as clut_jacobian, except calculations are based on the outer gamut of the CLUT. This function is used for linear extrapolation, which is enabled by the 'ubox' flag. The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMM function.

    GSA - a N-length vector (grid size array)

    X - a N-length input vector

    CACHE - a P-by-M matrix (CLUT) in cached format (where P is total number of grid points)

    J - an M x N Jacobian matrix

Example

    use ICC::Profile;
    use ICC::Support::Lapack;
    use Data::Dumper;

    $profile = ICC::Profile->new('../../test_data/profiles/originals/GRACoL2006_FullK_i1P.icc'); # open a profile
    $A2B0 = $profile->tag('A2B0');
    $clut = $A2B0->clut(); # get clut tag
    $cache = ICC::Support::Lapack::cache_2D($clut->array()); # get cached CLUT
    $gsa = $clut->gsa(); # get gsa

    $x = [0.1, 0.3, 0.5, 0.1]; # input vector
    $jac = ICC::Support::Lapack::clut_jacobian_ext($gsa, $x, $cache); # Jacobian matrix
    print Dumper($x, $jac);

Functions Used By matf.pm Module

The 'matf.pm' module implements a matrix transform with offset:

    Y = AX + B

In Perl, it is most efficient to store an array of vectors is like this:

    $X = [
        [vector 1],
        [vector 2],
           ...
        [vector K],
    ];

This matrix must be transposed before calculating Y = AX + B, and the result (Y) transposed as well. This transposing may be avoided by using the following identity:

    Yᵀ = XᵀAᵀ + Bᵀ

matf_fit

matf_fit determines optimum values for A and B from a dataset of input and output values, using least squares. Uses the LAPACK DGELS function.

Arguments

    X - a K-by-N matrix (input vectors)

    Y - a K-by-M matrix (output vectors)

    OFFSET - offset flag (use offset if true)

    INFO - return value (see DGELS)

    AB - solution containing M-by-N matrix A and, optionally, M-length vector B (see DGELS)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $x = [[1, 2, 3], [8, 3, 1], [5, 9, -1], [3, 3, 7], [7, 9, 2]];
    $y = [[3], [7], [11], [0], [8]];

    ($info, $ab) = ICC::Support::Lapack::matf_fit($x, $y, 0); # without offset
    print Dumper($info, $ab);

    ($info, $ab) = ICC::Support::Lapack::matf_fit($x, $y, -1); # with offset
    print Dumper($info, $ab);

matf_mat_trans

matf_mat_trans multiplies an array of vectors by a matrix, and adds an offset vector.

    Yᵀ = AXᵀ + B

Arguments

    A - a M-by-N matrix

    B - a M-length vector (offset), zeros if undefined or empty.

    X - a K-by-N matrix (input vectors)

    Y - a K-by-M matrix (output vectors)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1], [-1, 11, 12]];
    $b = [0.1, 0.2, 0.3, 0.4];
    $x = [[1, 2, 3], [0, 0, 0]];
    $y = ICC::Support::Lapack::matf_mat_trans($x, $a, $b);
    print Dumper($y);

matf_vec_trans

matf_vec_trans multiplies a single vector by a matrix, and adds an offset vector.

    Y = AX + B

Arguments

    A - a M-by-N matrix

    B - a M-length vector (offset), zeros if undefined or empty.

    X - a N-length vector (input)

    Y - a M-length vector (output)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1], [-1, 11, 12]];
    $b = [0.1, 0.2, 0.3, 0.4];
    $x = [1, 2, 3];
    $y = ICC::Support::Lapack::matf_vec_trans($x, $a, $b);
    print Dumper($y);

matf_inv

matf_inv solves the following equation for X, given Y, A and B. the matrix A must be square.

    Yᵀ = AXᵀ + B

Arguments

    A - a M-by-M matrix (square)

    B - a M-length vector (offset), zeros if undefined or empty.

    X - a K-by-M matrix (input vectors)

    Y - a K-by-M matrix (output vectors)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [4, -1, 6], [7, 8, -1]];
    $b = [0.1, 0.2, 0.3];
    $y = [[1, 2, 3], [0, 0, 0]];
    ($info, $x) = ICC::Support::Lapack::matf_inv($y, $a, $b);
    print Dumper($x);

Functions Used By nMIX.pm Module

The 'nMIX.pm' module implements a spectral ink mixing model based on the Yule-Nielsen equation. The model supports any number of input channels (inks) and any number of output channels (spectral values). The corner points of the process (device values equal zero or one) are stored in the object. These spectral values are exponentiated by a 'delta' value, mixed using barycentric coordinates computed from the device values, and exponentiated by the reciprocal delta value.

nMIX_mat_trans

nMIX_mat_trans computes a list of output vectors (matrix) from a list of input vectors (matrix). The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMM function.

Arguments

    X - a K-by-N matrix (input vectors)

    CACHE - a P-by-M matrix (coefficients) in cached format (where P is size of barycentric coordinate vector)

    D - an optional M-length exponent vector

    Y - a K-by-M matrix (output vectors)

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use ICC::Support::nMIX;
    use Data::Dumper;

    $dataset = '~/Documents/ProfileMaker/Measurements/FOGRA/2015_FOGRA51_FOGRA52/CGATS_format/FOGRA51_Spectral_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object
    $x = $chart->device([]); # get device values

    $nMIX = ICC::Support::nMIX->new($chart); # make nMIX object
    $array = $nMIX->array(); # get corner point array
    $cache = ICC::Support::Lapack::cache_2D($array->power(0.5)); # make cached version of corner point array

    $d = [(2) x @{$array->[0]}]; # make inverse delta vector (1/0.5 = 2)

    $x = [[0, 0, 0, 0], [0.5, 0.5, 0.5, 0], [0.1, 0.2, 0.3, 0.4]]; # matrix of input vectors
    $y = ICC::Support::Lapack::nMIX_mat_trans($x, $cache, $d); # transform input matrix
    print Dumper($x, $y);

nMIX_vec_trans

nMIX_vec_trans computes a single output vector from an input vector. The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMV function.

Arguments

    X - a N-length input vector

    CACHE - a P-by-M matrix (coefficients) in cached format (where P is size of barycentric coordinate vector)

    D - an optional M-length exponent vector

    Y - a M-length output vector

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use ICC::Support::nMIX;
    use Data::Dumper;

    $dataset = '~/Documents/ProfileMaker/Measurements/FOGRA/2015_FOGRA51_FOGRA52/CGATS_format/FOGRA51_Spectral_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object
    $x = $chart->device([]); # get device values

    $nMIX = ICC::Support::nMIX->new($chart); # make nMIX object
    $array = $nMIX->array(); # get corner point array
    $cache = ICC::Support::Lapack::cache_2D($array->power(0.5)); # make cached version of corner point array

    $d = [(2) x @{$array->[0]}]; # make inverse delta vector (1/0.5 = 2)

    $x = [0.1, 0.2, 0.3, 0.4]; # input vector
    $y = ICC::Support::Lapack::nMIX_vec_trans($x, $cache, $d); # transform input vector
    print Dumper($x, $y);

nMIX_jacobian

nMIX_jacobian computes the Jacobian matrix for a single sample from its input and output vectors. The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMM function.

Arguments

    X - a N-length input vector

    CACHE - a P-by-M matrix (coefficients) in cached format (where P is size of barycentric coordinate vector)

    D - an optional M-length exponent vector

    Y - an optional M-length output vector (required if D is supplied)

    J - an M x N Jacobian matrix

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use ICC::Support::nMIX;
    use Data::Dumper;

    $dataset = '~/Documents/ProfileMaker/Measurements/FOGRA/2015_FOGRA51_FOGRA52/CGATS_format/FOGRA51_Spectral_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object
    $x = $chart->device([]); # get device values

    $nMIX = ICC::Support::nMIX->new($chart); # make nMIX object
    $array = $nMIX->array(); # get corner point array
    $cache = ICC::Support::Lapack::cache_2D($array->power(0.5)); # make cached version of corner point array

    $d = [(2) x @{$array->[0]}]; # make inverse delta vector (1/0.5 = 2)

    $x = [0.1, 0.2, 0.3, 0.4]; # input vector
    $y = ICC::Support::Lapack::nMIX_vec_trans($x, $cache, $d); # transform input vector
    $jac = ICC::Support::Lapack::nMIX_jacobian($x, $cache, $d, $y); # compute Jacobian matrix
    print Dumper($x, $jac, $y);

nMIX_power

nMIX_power exponentiates a matrix to a vector of exponents. Matrix values are tested for sign.

Arguments

    X - an M x N input matrix

    D - an N-length exponent vector

    Y - an M x N output matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $x = [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]; 
    $d = [1, 2, 3, 4];

    $y = ICC::Support::Lapack::nMIX_power($x, $d);
    print Dumper($x, $y);

Functions Used By nNET.pm Module

The 'nNET.pm' module implements a simple feed-forward neural network. The neural network has a variable number of nodes which may be either radial basis functions, or Perl CODE blocks. Each node receives all of the inputs, and outputs a single value. These output values become the input to a matrix multiply function, with optional offset. The neural network is 'trained' by fitting a dataset of input and output values. The fitting data will often be rank deficient, which would cause LAPACK DGELS to fail. Instead, we call DGELSS, which uses the singular value decomposition (SVD). DGELSS fits rank deficient data reliably.

nNET_fit

nNET_fit determines optimum values for A and B from a dataset of input and output values, using least squares. Uses the LAPACK DGELSS function, which handles a rank-deficient X matrix.

Arguments

    X - a K-by-N matrix (input vectors)

    Y - a K-by-M matrix (output vectors)

    OFFSET - offset flag (use offset if true)

    INFO - return value (see DGELSS)

    AB - solution containing M-by-N matrix A and, optionally, M-length vector B (see DGELSS)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $x = [[1, 2, 3], [8, 3, 1], [5, 9, -1], [3, 3, 7], [7, 9, 2]];
    $y = [[3], [7], [11], [0], [8]];

    ($info, $ab) = ICC::Support::Lapack::nNET_fit($x, $y, 0); # without offset
    print Dumper($info, $ab);

    ($info, $ab) = ICC::Support::Lapack::nNET_fit($x, $y, -1); # with offset
    print Dumper($info, $ab);

Functions Used By nPINT.pm Module

The 'nPINT.pm' module implements a multi-dimensional polynomial interpolator. It supports any number of inputs and outputs, with variable polynomial degree for each input. The nPINT module is similar to the matf module, except that the input values are converted to an 'extended' form, based on the 'pda' array, which contains the polynomial degree of each input. This is explained in more detail by the document 'nPINT_module_notes.txt'.

nPINT_fit

nPINT_fit determines optimum coefficients from a dataset of input and output values, using least squares. Uses the LAPACK DGELS function.

Arguments

    PDA - a N-length vector (polynomial degree array)

    X - a K-by-N matrix (input vectors)

    Y - a K-by-M matrix (output vectors)

    INFO - return value (see DGELS)

    COEF - a X-by-M optimum coefficient matrix (where X is size of extended input vector)

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use Data::Dumper;

    $dataset = '../../test_data/misc/GRACoL2006_Coated1_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object from GRACoL2006 dataset
    $pda = [4, 4, 4, 4]; # polynomial degree for each channel
    $x = $chart->device([]); # get device values
    $y = $chart->lab([]); # get L*a*b* values
    ($info, $coef) = ICC::Support::Lapack::nPINT_fit($pda, $x, $y);
    print Dumper($info, $coef);

nPINT_mat_trans

nPINT_mat_trans computes a matrix of extended device values, and mutliplies it by the coefficient matrix, returning a matrix of output values. The coefficent matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMM function.

Arguments

    PDA - a N-length vector (polynomial degree array)

    X - a K-by-N matrix (input vectors)

    CACHE - a P-by-M matrix (coefficients) in cached format (where P is size of extended input vector)

    Y - a K-by-M matrix (output vectors)

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use Data::Dumper;

    $dataset = '~/Documents/ProfileMaker/Measurements/GRACoL/GRACoL2006_Release/GRACoL2006_Coated1_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object from GRACoL2006 dataset
    $pda = [4, 4, 4, 4]; # polynomial degree for each channel
    $x = $chart->device([]); # get device values
    $y = $chart->lab([]); # get L*a*b* values
    ($info, $coef) = ICC::Support::Lapack::nPINT_fit($pda, $x, $y);

    $cols = ICC::Support::Lapack::xcols($pda); # get number of coefficients
    $cache = ICC::Support::Lapack::cache_2D([@{$coef}[0 .. $cols - 1]]); # cache the coefficient matrix

    $x = [[0, 0, 0, 0], [0.5, 0.5, 0.5, 0], [0, 0, 0, 1]]; # matrix of input vectors
    $y = ICC::Support::Lapack::nPINT_mat_trans($pda, $x, $cache);
    print Dumper($x, $y);

nPINT_vec_trans

nPINT_vec_trans computes a vector of extended device values, and multiplies it by the coefficient matrix, returning an output vector. The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMV function.

Arguments

    PDA - a N-length vector (polynomial degree array)

    X - a N-length input vector

    CACHE - a P-by-M matrix (coefficients) in cached format (where P is size of extended input vector)

    Y - a M-length output vector

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use Data::Dumper;

    $dataset = '~/Documents/ProfileMaker/Measurements/GRACoL/GRACoL2006_Release/GRACoL2006_Coated1_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object from GRACoL2006 dataset
    $pda = [4, 4, 4, 4]; # polynomial degree for each channel
    $x = $chart->device([]); # get device values
    $y = $chart->lab([]); # get L*a*b* values
    ($info, $coef) = ICC::Support::Lapack::nPINT_fit($pda, $x, $y);

    $cols = ICC::Support::Lapack::xcols($pda); # get number of coefficients
    $cache = ICC::Support::Lapack::cache_2D([@{$coef}[0 .. $cols - 1]]); # cache the coefficient matrix

    $x = [0, 0, 0, 0]; # input vector
    $y = ICC::Support::Lapack::nPINT_vec_trans($pda, $x, $cache);
    print Dumper($x, $y);

nPINT_jacobian

nPINT_jacobian computes the Jacobian matrix from an input vector. The coefficient matrix is supplied in cached format, to improve speed. Uses the BLAS DGEMV function.

Arguments

    PDA - a N-length vector (polynomial degree array)

    X - a N-length input vector

    CACHE - a P-by-M matrix (coefficients) in cached format (where P is size of extended input vector)

    JAC - a M-by-N matrix (Jacobian matrix)

Example

    use ICC::Support::Lapack;
    use ICC::Support::Chart;
    use Data::Dumper;

    $dataset = '~/Documents/ProfileMaker/Measurements/GRACoL/GRACoL2006_Release/GRACoL2006_Coated1_CGATS.txt';
    $chart = ICC::Support::Chart->new($dataset); # make chart object from GRACoL2006 dataset
    $pda = [4, 4, 4, 4]; # polynomial degree for each channel
    $x = $chart->device([]); # get device values
    $y = $chart->lab([]); # get L*a*b* values
    ($info, $coef) = ICC::Support::Lapack::nPINT_fit($pda, $x, $y);

    $cols = ICC::Support::Lapack::xcols($pda); # get number of coefficients
    $cache = ICC::Support::Lapack::cache_2D([@{$coef}[0 .. $cols - 1]]); # cache the coefficient matrix

    $x = [0, 0, 0, 0]; # input vector
    $jac = ICC::Support::Lapack::nPINT_jacobian($pda, $x, $cache);
    print Dumper($x, $jac);

xcols

xcols computes the number of 'extended device values' from the PDA array. This number is the product of each array element plus one. For example, if the PDA array is [4, 4, 4, 4], the number is 5 x 5 x 5 x 5 = 625.

Arguments

    PDA - a N-length vector (polynomial degree array)

    COLS - number of extended device values

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $pda = [4, 4, 4, 4]; # polynomial degree for each channel
    $cols = ICC::Support::Lapack::xcols($pda); # get number of coefficients
    print Dumper($pda, $cols);

cache_2D

cache_2D converts a Perl 2-D array into a double C array, returned as a Perl string. The C array is in column-major order (FORTRAN) as used by LAPACK and BLAS. The cache is a reference to a Perl array containing this string and the number of rows and columns.

Arguments

    COEF - a numeric Perl 2-D array

    CACHE - a Perl array containing the cache string, the number of rows, and the number of columns

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $coef = [[0, 0, 0], [1, 2, 3], [4, 5, 6], [7, 8, 9]];
    $cache = ICC::Support::Lapack::cache_2D($coef);
    print Dumper($coef, $cache);

Matrix Functions

Some commonly used matrix functions are provided. These functions use the LAPACK library for improved speed and accuracy (compared to a pure PERL implementation).

zeros

Creates a matrix (2-D array) of 0's.

Arguments

    M - number of rows

    N - number of columns

    ZEROS - a M-by-N matrix of 0's

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $m = 4;
    $n = 3;

    # creates a square matrix with one parameter
    $zeros = ICC::Support::Lapack::zeros($m);
    print Dumper($zeros);

    # creates a rectangular matrix with two parameters
    $zeros = ICC::Support::Lapack::zeros($m, $n);
    print Dumper($zeros);

ones

Creates a matrix (2-D array) of 1's.

Arguments

    M - number of rows

    N - number of columns

    ONES - a M-by-N matrix of 1's

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $m = 4;
    $n = 3;

    # creates a square matrix with one parameter
    $ones = ICC::Support::Lapack::ones($m);
    print Dumper($ones);

    # creates a rectangular matrix with two parameters
    $ones = ICC::Support::Lapack::ones($m, $n);
    print Dumper($ones);

eye

Creates an identity matrix (2-D array).

Arguments

    M - number of rows

    N - number of columns

    IDN - a M-by-N identity matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $m = 4;
    $n = 3;

    # creates a square matrix with one parameter
    $idn = ICC::Support::Lapack::eye($m);
    print Dumper($idn);

    # creates a rectangular matrix with two parameters
    $idn = ICC::Support::Lapack::eye($m, $n);
    print Dumper($idn);

power

Exponentiates a matrix. The exponent may be a scalar or a vector. The vector size is the column size of the matrix.

Arguments

    X - a M-by-N matrix

    S - a numeric scalar

    V - a vector of N elements

    Y - a M-by-N matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $x = [[0, 1, 2], [3, 4, 5]]; # matrix

    $v = [1, 2, 3]; # exponent vector

    $y = ICC::Support::Lapack::power($x, 2); # square all matrix elements
    print Dumper($x, $y);

    $y = ICC::Support::Lapack::power($x, $v); # square second column, cube third column
    print Dumper($x, $y);

chol

Computes Cholesky factorization of a matrix using the LAPACK subroutine DPOTRF. The input matrix must be symmetric and positive-definite.

Arguments

    A - a N-by-N matrix

    INFO - return value (see DPOTRF)

    CHOL - a N-by-N matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[13, -6, 5], [-6, 8, -2], [5, -2, 5]];
    ($info, $chol) = ICC::Support::Lapack::chol($a);
    print Dumper($a, $info, $chol);

solve

Solves the matrix equation AX = B using the LAPACK subroutine DGESV.

Arguments

    A - a N-by-N matrix

    B - a N-by-NRHS matrix

    INFO - return value (see DGESV)

    X - a N-by-NRHS matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[3, 1, 3], [1, 5, 9], [2, 6, 5]];
    $b = [[1], [-2], [3]];
    ($info, $x) = ICC::Support::Lapack::solve($a, $b);
    print Dumper($a, $b, $info, $x);

trisolve

Solves the matrix equation AX = B using the LAPACK subroutine DGTSV. A is a tridiagonal matrix represented by three vectors, DL, D, and DU.

Arguments

    DL - a N-1 vector

    D - a N vector

    DU - a N-1 vector

    B - a N-by-NRHS matrix

    INFO - return value (see DGTSV)

    X - a N-by-NRHS matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $dl = [1, 1, 1, 1, 1]; # lower diagonal
    $d = [2, 4, 4, 4, 4, 2]; # diagonal
    $du = [1, 1, 1, 1, 1]; # upper diagonal
    $b = [[1], [-2], [3], [0], [-5], [8]];
    ($info, $x) = ICC::Support::Lapack::trisolve($dl, $d, $du, $b);
    print Dumper($dl, $d, $du, $b, $info, $x);

inv

Inverts a square matrix. Uses LAPACK subroutines DGETRF and DGETRI. If an error occurs, the name of the subroutine that failed is returned in X.

Arguments

    A - a N-by-N matrix

    INFO - return value (see DGETRF and DGETRI)

    X - a N-by-N matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[3, 1, 3], [1, 5, 9], [2, 6, 5]];
    ($info, $x) = ICC::Support::Lapack::inv($a);
    print Dumper($a, $info, $x);

normal

Solves the normal equations using QR factorization, with optional sample weighting. Uses LAPACK subroutines DGEQRF, DORMQR and DTRTRS. If an error occurs, the name of the subroutine that failed is returned in DELTA.

Arguments

    A - a K-by-N Jacobian matrix

    B - a K-by-1 residual matrix

    W - a K-by-1 weight matrix (optional)

    INFO - return value (see DGEQRF, DORMQR and DTRTRS)

    DELTA - a N-by-1 delta matrix (or the name of failed subroutine)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[3, 1, 3], [1, 5, 9], [2, 6, 5], [7, -2, 8]];
    $b = [[1], [-2], [3], [-1]];
    $w = [[1], [1], [0.5], [0.5]];

    ($info, $delta) = ICC::Support::Lapack::normal($a, $b); # uniform weights
    print Dumper($a, $b, $info, $delta);

    ($info, $delta) = ICC::Support::Lapack::normal($a, $b, $w); # optional weight matrix
    print Dumper($a, $b, $w, $info, $delta);

qr

Computes a QR factorization of a matrix. Uses LAPACK subroutines DGEQRF and DORMQR. If an error occurs, the name of the subroutine that failed is returned in Q, and R is undefined. The economy flag, if set to zero, reduces the intermediate dimension when m > n.

Arguments

    A - a M-by-N input matrix

    E - economy flag, trims output when set to zero

    INFO - return value (see DGEQRF and DORMQR)

    Q - a M by M -or- M by N orthogonal matrix

    R - a M-by-N -or- N by N right triangular matrix

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [4, 11, 6], [7, 8, 9]]; # m = n
    ($info, $q, $r) = ICC::Support::Lapack::qr($a);
    print Dumper($a, $info, $q, $r);

    $a = [[1, 2, 3], [4, 11, 6], [7, 8, 9], [10, 4, 11]]; # m > n
    ($info, $q, $r) = ICC::Support::Lapack::qr($a);
    print Dumper($a, $info, $q, $r);

    $a = [[1, 2, 3, 4], [5, 6, 11, 8], [8, 6, 10, 11]];; # m < n
    ($info, $q, $r) = ICC::Support::Lapack::qr($a);
    print Dumper($a, $info, $q, $r);

    $a = [[1, 2, 3], [4, 11, 6], [7, 8, 9], [10, 4, 11]]; # m > n
    ($info, $q, $r) = ICC::Support::Lapack::qr($a, 0); # economy flag
    print Dumper($a, $info, $q, $r);

dot

Computes the dot product of two vectors using the BLAS subroutine DDOT.

    P = V1 • V2

Arguments

    V1 - a vector (dimension N)

    V2 - a vector (dimension N)

    P - a scalar

Example

    use ICC::Support::Lapack;

    $a = [1, 2, 3, 4, 5, 6];
    $b = [6, 5, 4, 3, 2, 1];

    $prod = ICC::Support::Lapack::dot($a, $b);

    print $prod;

vec_xplus

Multiplies a matrix and a vector, adding an optional offset vector, using the BLAS subroutine DGEMV.

    Y = A * X + C

The DGEMV parameters TRANS, ALPHA and BETA may be set with an optional hash parameter. Default value for TRANS is 'N'. Default values for ALPHA and BETA are 1.0.

Arguments

    A - an M-by-N matrix

    X - a vector, dimension N or M (depends on TRANS)

    C - a vector, optional offset, dimension M or N (depends on TRANS)

    HASH - parameter hash reference (optional)

    Y - a vector, result, dimension M or N (depends on TRANS)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 7, -2, 4], [0, 5, -3, 1], [-4, 8, 9, 3]];
    $at = [[1, 0, -4], [7, 5, 8], [-2, -3, 9], [4, 1, 3]];
    $v1 = [6, 5, 4, 3];
    $v2 = [1, 2, 3];
    $c1 = [0.1, 0.2, 0.3];
    $c2 = [0.1, 0.2, 0.3, 0.4];

    $y = ICC::Support::Lapack::vec_xplus($a, $v1);
    print Dumper($y);

    $y = ICC::Support::Lapack::vec_xplus($a, $v1, $c1); # with offset
    print Dumper($y);

    $y = ICC::Support::Lapack::vec_xplus($a, $v1, {'alpha' => 2}); # multiply product by 2
    print Dumper($y);

    $y = ICC::Support::Lapack::vec_xplus($a, $v1, $c1, {'beta' => 2}); # multiply offset by 2
    print Dumper($y);

    $y = ICC::Support::Lapack::vec_xplus($at, $v1, {'trans' => 'T'}); # transpose
    print Dumper($y);

    $y = ICC::Support::Lapack::vec_xplus($a, $v2, {'trans' => 'T'}); # transpose
    print Dumper($y);

    $y = ICC::Support::Lapack::vec_xplus($a, $v2, $c2, {'trans' => 'T'}); # transpose with offset
    print Dumper($y);

mat_xplus

Multiplies two matrices and adds optional offset matrix, using the BLAS subroutine DGEMM.

    Y = A * B + C

The DGEMM parameters TRANSA, TRANSB, ALPHA and BETA may be set with an optional hash parameter. Default values for TRANSA and TRANSB are 'N'. Default values for ALPHA and BETA are 1.0.

Arguments

    A - an M-by-K or K-by-M matrix (depends on TRANSA)

    B - a K-by-N or N-by-K matrix (depends on TRANSB)

    C - an M-by-N matrix (optional)

    HASH - parameter hash reference (optional)

    Y - an M-by-N matrix (result)

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [8, 3, 1], [5, 9, -1], [3, 3, 7], [7, 9, 2]];
    $at = [[1, 8, 5, 3, 7], [2, 3, 9, 3, 9], [3, 1, -1, 7, 2]]; # transposed
    $b = [[3, 5], [7, -2], [11, 0]];
    $bt = [[3, 7, 11], [5, -2, 0]]; # transposed
    $c = [[0, 0.1], [0.2, 0.3], [0.4, 0.5], [0.6, 0.7], [0.8, 0.9]];

    $y = ICC::Support::Lapack::mat_xplus($a, $b); # no offset
    print Dumper($y);

    $y = ICC::Support::Lapack::mat_xplus($a, $b, $c); # with offset
    print Dumper($y);

    $y = ICC::Support::Lapack::mat_xplus($at, $bt, {'transa' => 'T', 'transb' => 'T'}); # transpose
    print Dumper($y);

    $y = ICC::Support::Lapack::mat_xplus($at, $bt, $c, {'transa' => 'T', 'transb' => 'T'}); # transpose with offset
    print Dumper($y);

    $y = ICC::Support::Lapack::mat_xplus($b, $a, {'transa' => 'T', 'transb' => 'T'}); # transpose and swap $a and $b
    print Dumper($y);

    $y = ICC::Support::Lapack::mat_xplus($a, $b, {'alpha' => 2}); # multiply product by 2
    print Dumper($y);

    $y = ICC::Support::Lapack::mat_xplus($a, $b, $c, {'beta' => 2}); # multiply offset by 2
    print Dumper($y);

Statistics Functions

Some commonly used statistics functions are provided. These functions use the LAPACK library for improved speed and accuracy (compared to a pure PERL implementation).

mean

Computes the mean values of a sample set.

Arguments

    A - a M-by-N sample matrix

    MEAN - a N-element array containing the column mean values

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 3], [7, 5, 6], [10, 8, 9], [4, 11, 12]];
    $mean = ICC::Support::Lapack::mean($a);
    print Dumper($a, $mean);

cov

Computes the covariance matrix of a sample set. Uses the BLAS subroutine DGEMM

Arguments

    A - a M-by-N sample matrix

    COV - a N-by-N covariance matrix

    MEAN - a N-element array containing the column mean values

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 9], [7, 5, 12], [10, 8, 3], [4, 11, 6]];
    ($cov, $mean) = ICC::Support::Lapack::cov($a);
    print Dumper($a, $cov, $mean);

pca

Computes a principal component analysis of a sample set. Uses the LAPACK subroutine DGESVD.

Arguments

    A - a M-by-N sample matrix

    INFO - return value (see DGESVD)

    U - a M-by-M orthogonal matrix

    SIGMA - a M-by-N diagonal matrix (an array of the diagonal values)

    VT - a N-by-N orthogonal matrix (transposed)

    MEAN - a N-element array containing the column mean values

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $a = [[1, 2, 9], [7, 5, 12], [10, 8, 3], [4, 11, 6]];
    ($info, $u, $sigma, $vt, $mean) = ICC::Support::Lapack::pca($a);
    print Dumper($a, $info, $u, $sigma, $vt, $mean);

mahal

Computes the Mahalanobis distance between two vectors.

Arguments

    All values are real

    V1 - a N-length vector

    V2 - a N-length vector

    COV - a N-by-N covariance matrix

    COVINV - a N-by-N inverse of covariance matrix

    DIST - the Mahalanobis distance between the vectors

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $v1 = [1, 2, 3];
    $v2 = [6, 5, 4];
    $cov = [[13, -6, 5], [-6, 8, -2], [5, -2, 5]];

    ($info, $covinv) = ICC::Support::Lapack::inv($cov);
    ($info == 0) || die('inverting covariance failed');

    $dist = ICC::Support::Lapack::mahal($v1, $v2, $covinv);
    print Dumper($v1, $v2, $cov, $covinv, $dist);

euclid

Computes the Euclidean distance between two vectors.

Arguments

    All values are real

    V1 - a N-length vector

    V2 - a N-length vector

    DIST - the Euclidean distance between the vectors

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $v1 = [1, 2, 3];
    $v2 = [6, 5, 4];

    $dist = ICC::Support::Lapack::euclid($v1, $v2);
    print Dumper($v1, $v2, $dist);

skew_normal

Computes the skew-normal probability distribution function (pdf). See http://azzalini.stat.unipd.it/SN/Intro/intro.html

Arguments

    All values are real

    X - input value

    LOCATION - input offset, ξ

    SCALE - steepness, ω, similar to sigma

    SHAPE - shape, α, if > 0, pdf is skewed right, if < 0, pdf is skewed left

    SN - output value

Example

    use ICC::Support::Lapack;
    use Data::Dumper;

    $x = 0;
    $location = 0;
    $scale = 1;
    $shape = 0.5;

    # compute skew_normal pdf
    $sn = ICC::Support::Lapack::skew_normal($x, $location, $scale, $shape);

    print "$x $sn\n";

Special Functions

expm1

Computes the expm1 function using the math.c library. This function is exp(arg) - 1.0, but more accurate for arguments near zero.

Arguments

    X - a real value

    EXP - a real value

Example

    use ICC::Support::Lapack;

    $x = 1.0;
    $exp = ICC::Support::Lapack::expm1($x);
    print "$x $exp\n";

log1p

Computes the log1p function using the math.c library. This function is log(1.0 + arg), but more accurate for arguments near zero.

Arguments

    X - a real value

    LOG - a real value

Example

    use ICC::Support::Lapack;

    $x = 1.0;
    $log = ICC::Support::Lapack::log1p($x);
    print "$x $log\n";

erf

Computes the error function using the math.c library.

Arguments

    X - a real value

    ERF - a real value

Example

    use ICC::Support::Lapack;

    $x = 1.0;
    $erf = ICC::Support::Lapack::erf($x);
    print "$x $erf\n";

tgamma

Computes the gamma function using the math.c library.

Arguments

    X - a real value

    GAMMA - a real value

Example

    use ICC::Support::Lapack;

    $x = 1.0;
    $gamma = ICC::Support::Lapack::tgamma($x);
    print "$x $gamma\n";

beta

Computes the beta function using the math.c library.

Arguments

    X - a real value
    Y - a real value

    BETA - a real value

Example

    use ICC::Support::Lapack;

    $x = 1.0;
    $y = 3.5
    $beta = ICC::Support::Lapack::beta($x, $y);
    print "$x $y $beta\n";

Miscellaneous Functions

round

Computes the round function using the math.c library.

Arguments

    X - a real value

    ROUND - a real value

Example

    use ICC::Support::Lapack;

    $x = 2.4;
    $round = ICC::Support::Lapack::round($x);
    print "$x $round\n";

bandpass

Computes a bandpass function, y = (1/sigma) * exp(-(abs(x - mu)/sigma)^gamma). The function is optionally asymmetric, using sigma1 and gamma1 when x < mu, sigma2 and gamma2 when x >= mu. Parameters sigma2, gamma1 and gamma2 are optional. See examples below for parameter defaults. The total integral of the function is equal to 1. Sigma values must not be zero. Gamma values must be > 1/171. With three parameters, the bandpass function is the same as the normal (Gaussian) pdf.

Arguments

    All values are real

    X - input value

    MU - offset
    SIG1 - sigma (X < MU)
    SIG2 - sigma (X >= MU)
    GAM1 - gamma (X < MU)
    GAM2 - gamma (X >= MU)

    Y - output value

Example

    use ICC::Support::Lapack;

    $mu = 0.5;
    $sig1 = 1;
    $sig2 = 1.5;
    $gam1 = 2;
    $gam2 = 2.3;

    $x = 0.25; # input value

    $y = ICC::Support::Lapack::bandpass($x, $mu, $sig1); # sig2 = sig1, gam1 = gam2 = 2
    print "$x $y\n";

    $y = ICC::Support::Lapack::bandpass($x, $mu, $sig1, $sig2); # gam1 = gam2 = 2
    print "$x $y\n";

    $y = ICC::Support::Lapack::bandpass($x, $mu, $sig1, $sig2, $gam1); # gam2 = gam1
    print "$x $y\n";

    $y = ICC::Support::Lapack::bandpass($x, $mu, $sig1, $sig2, $gam1, $gam2);
    print "$x $y\n";

hermite

Computes cubic spline value using Hermite form.

Arguments

    All values are real

    T - interval position (0 - 1)

    Y0 - value at position 0
    Y1 - value at position 1
    M0 - derivative at position 0
    M1 - derivative at position 1

    Y - interpolated value

Example

    use ICC::Support::Lapack;

    $t = 0.25; # position in unit interval

    $y0 = 1; # y(0)
    $y1 = 2; # y(1)
    $m0 = 1.5; # y'(0)
    $m1 = 0.5; # y'(1)

    $y = ICC::Support::Lapack::hermite($t, $y0, $y1, $m0, $m1);
    print "$t $y\n";

bernstein

Computes Bernstein polynomial value. The Bernstein degree is determined by the size ($#) of the coefficient array. The maximum Bernstein degree is 20.

Arguments

    All values are real

    T - interval position (0 - 1)

    COEF - Bernstein polynomial coefficients

    Y - interpolated value

Example

    use ICC::Support::Lapack;

    $t = 0.25; # position in unit interval

    $coef = [1, 2, 3]; # Bernstein coefficients

    $y = ICC::Support::Lapack::bernstein($t, $coef);
    print "$t $y\n";

LICENSE

Programs in this distribution, authored by William B. Birkett, are licensed under the GNU General Public License, v3.

See http://www.gnu.org/licenses/gpl.html for license details.

AUTHOR

William B. Birkett, <wbirkett@doplganger.com>

NAME

SYNOPSIS

DESCRIPTION

FUNCTIONS

Interface to Selected LAPACK Routines

dgels

dgesv

dgesvd

dgeqrf

Interface to Selected BLAS Routines

dgemm

dgemv

Functions Used By clut.pm Module

clut_vec_trans

clut_jacobian

clut_jacobian_ext

Functions Used By matf.pm Module

matf_fit

matf_mat_trans

matf_vec_trans

matf_inv

Functions Used By nMIX.pm Module

nMIX_mat_trans

nMIX_vec_trans

nMIX_jacobian

nMIX_power

Functions Used By nNET.pm Module

nNET_fit

Functions Used By nPINT.pm Module

nPINT_fit

nPINT_mat_trans

nPINT_vec_trans

nPINT_jacobian

xcols

cache_2D

Matrix Functions

zeros

ones

eye

power

chol

solve

trisolve

inv

normal

qr

dot

vec_xplus

mat_xplus

Statistics Functions

mean

cov

pca

mahal

euclid

skew_normal

Special Functions

expm1

log1p

erf

tgamma

beta

Miscellaneous Functions

round

bandpass

hermite

bernstein

SEE ALSO

LICENSE

AUTHOR

COPYRIGHT