la_impl.h

/* -*- mode: C; c-basic-offset: 4 -*- */
/* ex: set shiftwidth=4 tabstop=4 expandtab: */
/*
 * Copyright (c) 2011-2012, Georgia Tech Research Corporation
 * All rights reserved.
 *
 * Author(s): Neil T. Dantam <ntd@gatech.edu>
 * Georgia Tech Humanoid Robotics Lab
 * Under Direction of Prof. Mike Stilman <mstilman@cc.gatech.edu>
 *
 *
 * This file is provided under the following "BSD-style" License:
 *
 *
 *   Redistribution and use in source and binary forms, with or
 *   without modification, are permitted provided that the following
 *   conditions are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
 *   CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 *   INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 *   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 *   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
 *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 *   USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 *   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *   POSSIBILITY OF SUCH DAMAGE.
 *
 */


#include "amino/def.h"

AA_API void AA_NAME(la,transpose)
( size_t m, size_t n,
  const AA_TYPE *A, size_t lda,
  AA_TYPE *B, size_t ldb );

AA_FDEC(AA_TYPE, la, ssd,
        size_t n,
        const AA_TYPE *x, size_t incx,
        const AA_TYPE *y, size_t incy )


AA_API void AA_NAME(la,lerp)
( size_t n, AA_TYPE u,
  const AA_TYPE *v1, size_t inc1,
  const AA_TYPE *v2, size_t inc2,
  AA_TYPE *vu, size_t incu );



AA_API void AA_NAME(la,3spline_param)
( size_t n, AA_TYPE tf,
  const AA_TYPE *x1, size_t incx1,
  const AA_TYPE *dx1, size_t incdx1,
  const AA_TYPE *x2, size_t inc2,
  const AA_TYPE *dx2, size_t incdx2,
  AA_TYPE *a2, AA_TYPE *a3 );

AA_API void AA_NAME(la,3spline)
( size_t n, AA_TYPE tf,
  const AA_TYPE *x1, size_t incx1,
  const AA_TYPE *dx1, size_t incdx1,
  const AA_TYPE *a2, const AA_TYPE *a3,
  AA_TYPE *x, size_t incx,
  AA_TYPE *dx, size_t incdx,
  AA_TYPE *ddx, size_t incddx );

AA_API void AA_NAME(la,5spline_param)
( size_t n, AA_TYPE tf,
  const AA_TYPE *x1, size_t incx1,
  const AA_TYPE *dx1, size_t incdx1,
  const AA_TYPE *ddx1, size_t incddx1,
  const AA_TYPE *x2, size_t inc2,
  const AA_TYPE *dx2, size_t incdx2,
  const AA_TYPE *ddx2, size_t incddx2,
  AA_TYPE *a3, AA_TYPE *a4, AA_TYPE *a5 );

AA_API void AA_NAME(la,5spline)
( size_t n, AA_TYPE tf,
  const AA_TYPE *x1, size_t incx1,
  const AA_TYPE *dx1, size_t incdx1,
  const AA_TYPE *ddx1, size_t incddx1,
  const AA_TYPE *a3, const AA_TYPE *a4, const AA_TYPE *a5,
  AA_TYPE *x, size_t incx,
  AA_TYPE *dx, size_t incdx,
  AA_TYPE *ddx, size_t incddx );

AA_FDEC(AA_TYPE, la, vecstd,
        size_t n,
        const AA_TYPE *x, size_t incx,
        AA_TYPE mu);

AA_FDEC(void, la, colmean,
        size_t m, size_t n,
        const AA_TYPE *A, size_t lda,
        AA_TYPE *x);

AA_FDEC(void, la, rowmean,
        size_t m, size_t n,
        const AA_TYPE *A, size_t lda,
        AA_TYPE *x);

/* AA_API void AA_NAME(la,_cmean) */
/* ( size_t m, size_t n, */
/*   const AA_TYPE *A, size_t lda, */
/*   AA_TYPE *x ); */

AA_FDEC(void, la, colcov,
        size_t m, size_t n,
        const AA_TYPE *A, size_t lda,
        const AA_TYPE *x,
        AA_TYPE *E, size_t lde);

/* AA_API void AA_NAME(la,_ccov) */
/* ( size_t m, size_t n, */
/*   const AA_TYPE *A, size_t lda, */
/*   const AA_TYPE *x, */
/*   AA_TYPE *E, size_t lde ); */



/* Hungarian algorithm to solve min assignment problem
 * \param n rows and cols of A
 * \param A cost matrix, column major, destroyed on exit
 * \param[out] row_assign array of column assigned to row at index i,
 *             length m. Unmatched elements have the value of -1.
 * \param col_assign mapping from col to matched row
 * \param lda leading dimension of A
 * \param row_assign array of column assigned to row at index i
 * \param iwork array of size 3*n*n +2*n
 */
/* AA_API void AA_NAME(la,opt_hungarian) */
/* ( size_t n, AA_TYPE *A, size_t lda, */
/*   ssize_t *row_assign, */
/*   ssize_t *col_assign, */
/*   ssize_t *iwork); */


/* Compute minimum iwork size for hungarian algorithm
 *
 * \sa aa_la_d_opt_hungarian
 */
/* static inline size_t AA_NAME(la,opt_hungarian_iworksize) */
/* ( size_t n ) { */
/*     return 3*n*n + 2*n; */
/* } */

/* Hungarian algorithm for rectangular distance matrix by padding with zeros.
 *
 * \param m rows of A
 * \param n cols of A
 * \param A distance matrix for minimization problem
 * \param lda leading dimension of A
 * \param[out] row_assign array of column assigned to row at index i,
 *             length m. Unmatched elements have the value of -1.
 * \param[out] col_assign mapping from col to matched row
 * \param work work array of length max(m,n)**2
 * \param iwork integer work array of length (3*max(m,n)**2 + 4*max(m,n))
 */
/* AA_API void AA_NAME(la,opt_hungarian_pad) */
/* ( size_t m, size_t n, const AA_TYPE *A, size_t lda, */
/*   ssize_t *row_assign, */
/*   ssize_t *col_assign, */
/*   AA_TYPE *work, ssize_t *iwork); */

/* Compute minimum iwork size for padded hungarian algorithm
 *
 * \sa aa_la_d_opt_hungarian_pad
 */
/* static inline size_t AA_NAME(la,opt_hungarian_pad_iworksize) */
/* (size_t m,size_t n) { */
/*     size_t p = AA_MAX(m,n); */
/*     return AA_NAME(la,opt_hungarian_iworksize)(p) + 2*p; */
/* } */


/* Compute minimum work size for padded hungarian algorithm
 *
 * \sa aa_la_d_opt_hungarian_pad
 */
/* static inline size_t AA_NAME(la,opt_hungarian_pad_worksize) */
/* (size_t m,size_t n) { */
/*     size_t p = AA_MAX(m,n); */
/*     return p*p; */
/* } */


/* Converts max assignment to min assignment for Hungarian algorithm.
 * \param m rows of A
 * \param n cols of A
 * \param A cost matrix for max problem,
 *        converted to matrix for min proble on exit
 * \param lda leading dimension of A
 */
/* AA_API void AA_NAME(la,opt_hungarian_max2min) */
/* ( size_t m, size_t n, AA_TYPE *A, size_t lda ); */


AA_FDEC( void, la, assign_hungarian,
         size_t m, size_t n,
         const AA_TYPE *A, size_t lda,
         ssize_t *row_assign, ssize_t *col_assign );

AA_FDEC( void, la, assign_hungarian_max2min,
         size_t m, size_t n,
         AA_TYPE *A, size_t lda );


static inline size_t AA_NAME(la,minloc)
( size_t n, AA_TYPE *x, size_t incx ) {
    size_t imin = 0;
    AA_TYPE xmin = *x;
    for( size_t i = 1; i < n; i ++ ) {
        if( x[i*incx] < xmin ) {
            imin = i;
            xmin = x[i*incx];
        }
    }
    return imin;
}

static inline AA_TYPE AA_NAME(la,mat_max)
( size_t m, size_t n, AA_TYPE *A, size_t lda,
  size_t *pi, size_t *pj ) {
    size_t im = 0;
    size_t jm = 0;
    AA_TYPE xm = *A;
    for( size_t i = 0; i < m; i++) {
        for( size_t j = 0; j < n; j++) {
            if( AA_MATREF(A, lda, i, j ) > xm ) {
                im = i;
                jm = j;
                xm =  AA_MATREF(A, lda, i, j );
            }
        }
    }
    if( pi ) *pi = im;
    if( pj ) *pj = jm;
    return xm;
}


static inline size_t AA_NAME(la,maxloc)
( size_t n, AA_TYPE *x, size_t incx ) {
    size_t imax = 0;
    AA_TYPE xmax = *x;
    for( size_t i = 1; i < n; i ++ ) {
        if( x[i*incx] > xmax ) {
            imax = i;
            xmax = x[i*incx];
        }
    }
    return imax;
}

AA_FDEC(AA_TYPE, la, angle,
        size_t n, const AA_TYPE *x, size_t incx,
        const AA_TYPE *y, size_t incy)




AA_API int AA_NAME(la,lls)
( size_t m, size_t n, size_t p,
  const AA_TYPE *A, size_t lda,
  const AA_TYPE *b, size_t ldb,
  AA_TYPE *x, size_t ldx );


AA_API int AA_NAME(la,svd)
( size_t m, size_t n, const AA_TYPE *A, size_t lda,
  AA_TYPE *U, size_t ldu,
  AA_TYPE *S,
  AA_TYPE *Vt, size_t ldvt );


AA_API AA_TYPE AA_NAME(la,median)
( size_t n, const AA_TYPE *x, size_t incx );

AA_API AA_TYPE AA_NAME(la,nmedian)
( size_t n, AA_TYPE *x );


static inline AA_TYPE AA_NAME(la,nmedian_pop)
( size_t n, AA_TYPE *x )
{
    AA_TYPE u = AA_NAME(la,nmedian)(n, x);
    aa_mem_region_local_pop(x);
    return u;
}

AA_TYPE AA_NAME(la,mad)
( size_t n, const AA_TYPE u, const AA_TYPE *x, size_t incx );

AA_TYPE AA_NAME(la,mad2)
( size_t m, size_t n, const AA_TYPE *u, const AA_TYPE *A, size_t lda );

AA_API int AA_NAME(la,eev)
( size_t n, const AA_TYPE *A, size_t lda,
  AA_TYPE *wr,
  AA_TYPE *wi,
  AA_TYPE *Vl, size_t ldvl,
  AA_TYPE *Vr,  size_t ldvr );

AA_FDEC( void, la, colfit,
         size_t m, size_t n,
         const AA_TYPE *A, size_t lda, AA_TYPE *x );


/* Ordering / Sort comparison function */
int AA_NAME(la,compar)( const void *a, const void *b );

int AA_NAME(la, ls_duqu_cmat)
(const double *sigma, const double *theta, double *cm);


#include "amino/undef.h"