matrix.h

Go to the documentation of this file.

/* matrixlib/matrix.h.  Generated by configure.  */
/*
 * matrix.h
 *
 * Programmer(s): Chris Harvey, Michael Paluszek, David Wilson, Lavanya Swetharanyan, Wil Turner
 *
 * Copyright 1996-2005, Princeton Satellite Systems. All rights reserved.
 *
 * MatrixLib is a complete solution for implementing one and two dimensional matrices, to solve embedded control
 * and simulation problems in the C++ environment. Virtually any engineering application will require the use of
 * matrices - this library implements all of the matrix operations supported by MATLAB in an easy to use and highly
 * efficient C++ API. It greatly reduces the time and money needed to port algorithms developed in MATLAB to C++.
 *
 */

#ifndef __MATRIX__
#define __MATRIX__

#include <stdio.h>
#include "int_array.h"
#include <string.h>
#include <stdlib.h>
#include <math.h>

extern "C++"
{

typedef enum {
    ml_no_error = 0,                
    ml_bad_index = 1,               
    ml_size_mismatch = 2,           
    ml_non_invertible_matrix = 3,   
    ml_bad_file_pointer = 4,        
    ml_bad_range = 5,               
    ml_bad_operand = 6,             
    ml_bad_string_format = 7,       
    ml_invalid_buffer = 8,          
    ml_lapack_error = 9             
} ml_error_codes;

#define ML_FP_TOL 1e-08
#define ML_ARR_INCR 4

/* FUNCTIONS USING BLAS
 *
 * ml_matrix::operator*=
 * ml_matrix::dot()
 * ml_matrix::vmag()
 */

/* FUNCTIONS USING LAPACK
 *
 * ml_matrix::inv()
 * solve_ax_eq_b()
 * svd()
 */

/* Use the BLAS if detected. */
/* #undef ML_USE_BLAS */

/* Use the LAPACK if detected. */
/* #undef ML_USE_LAPACK */

/* Linux OS */
/* #undef PSS_OS_LINUX */

/* Mac OS X */
#define PSS_OS_MACOSX 

/* UNIX OS */
#define PSS_OS_UNIX 

#define ML_COLUMN_MAJOR 

#ifdef ML_COLUMN_MAJOR
    #define ml_matloc(mat, x, y, rrows, rcols) (mat + (y) * (rrows) + (x))
    #define ml_incr_row(ptr, rrows, rcols) (ptr++)
    #define ml_incr_col(ptr, rrows, rcols) (ptr+=rrows)
    #define ml_major_stride(rrows,rcols) (rrows)
    #define ml_minor_stride(rrows,rcols) (rcols)
    #define ml_major_length(nrows,ncols) (nrows)
    #define ml_minor_length(nrows,ncols) (ncols)
#else
    #define ml_matloc(mat, x, y, rrows, rcols) (mat + (x) * (rcols) + (y))
    #define ml_incr_row(ptr, rrows, rcols) (ptr+=rcols)
    #define ml_incr_col(ptr, rrows, rcols) (ptr++)
    #define ml_major_stride(rrows,rcols) (rcols)
    #define ml_minor_stride(rrows,rcols) (rrows)
    #define ml_major_length(nrows,ncols) (ncols)
    #define ml_minor_length(nrows,ncols) (nrows)
#endif

class ml_matrix
{
    public:
        #if 0
        #pragma mark Construction and I/O
        #endif


        ml_matrix();

        inline ml_matrix(const ml_matrix &m)
        {
            // straight memcpy of parameter m into the contents of "this"
            memcpy(this, &m, sizeof(ml_matrix));
            double *t = mat;
            if (using_base)
            {
                mat = base_mat;
            }
            else
            {
                // need a new mat, since memcpy doesn't copy pointers
                mat = new double[real_cols * real_rows];
                // memcpy the contents of m.mat to this.mat
                memcpy(mat, t, len*sizeof(double));
            }
            buff = NULL;
            buffsize = 0;
            seed = 0;
        }

        inline ml_matrix(int rows, int cols = 1)
        {
            // if parameters are less than 0, set them to 1, otherwise use what was passed in
            num_rows = (rows < 0) ? 1 : rows;
            num_cols = (cols < 0) ? 1 : cols;
            int copysize;
            if ((num_rows <= ML_ARR_INCR) && (num_cols <= ML_ARR_INCR))
            {
                real_cols = real_rows = ML_ARR_INCR;
                len = ML_ARR_INCR*ML_ARR_INCR;
                copysize = sizeof(double) * len;
                mat = base_mat;
                using_base = 1;
            }
            else
            {
                // real_rows and real_cols will always be the nearest multiple of 4 above the requested
                // number of rows and cols. this is to allow extra room to prevent reallocation on small resizes
                // of the matrix.
                real_rows = ((num_rows>>2)+1) << 2;
                real_cols = ((num_cols>>2)+1) << 2;
                // len is the entire buffersize of mat, including the padding for potential resizes
                len = real_rows * real_cols;
                copysize = sizeof(double) * len;
                // allocate a new double array for mat, the size of len
                mat = new double[len];
                using_base = 0;
            }
            // set all elements of mat to 0
            memset(mat, 0, copysize);
            buff = NULL;
            buffsize = 0;
            seed = 0;
            error = 0;
        }
        ml_matrix(int rows, int cols, double value);
        ml_matrix(int rows, int cols, double *v);
        ml_matrix(const char *s);
        ml_matrix(void *bf);
        ml_matrix(FILE *f, int *linenum = NULL);
        inline ~ml_matrix()
        {
            if (mat && !using_base)
            {
                delete [] mat;
            }
            if (buff)
            {
                free(buff);
            }
        }
        char *to_string(int human_readable = 0);
        void display(char *s = NULL);
        void *to_binary(unsigned int *size = NULL) const;

        #if 0
        #pragma mark Matrix Manipulation
        #endif


        double &operator()(int row, int col);

        inline double get(int row, int col) const
        {
            if (row < 1 || col < 1)
            {
                return mat[0];
            }
            if((row > num_rows) || (col > num_cols))
            {
                return mat[0];
            }
            //printf("getting %d, %d; rrows is %d, rcols is %d\n",row,col,real_rows,real_cols);
            return *ml_matloc(mat,row-1,col-1,real_rows,real_cols);
        }
        void resize(int rows, int cols);
        ml_matrix &append_identity();
        ml_matrix sub_matrix(int from_row, int to_row, int from_col, int to_col) const;
        ml_matrix sub_matrix(int n_rows, int *rows, int n_cols, int *cols) const;
        ml_matrix sub_matrix(const ml_int_array &rows, const ml_int_array &cols) const;
        ml_matrix get_row(int row) const;
        ml_matrix get_col(int col) const;
        ml_matrix get_row_set(int n_rows, int *rows) const;
        ml_matrix get_row_set(const ml_int_array &rows) const;
        ml_matrix get_col_set(int n_cols, int *cols) const;
        ml_matrix get_col_set(const ml_int_array &cols) const;
        ml_matrix sub_mat_row(int row, int n_cols, int *cols) const;
        ml_matrix sub_mat_row(int row, const ml_int_array &cols) const;
        ml_matrix sub_mat_col(int n_rows, int *rows, int col) const;
        ml_matrix sub_mat_col(const ml_int_array &rows, int col) const;
        ml_matrix extract_elements(const ml_matrix &elements);
        ml_matrix &delete_rows_cols(int n_rows, int *rows, int n_cols, int *cols);
        ml_matrix &delete_rows_cols(const ml_int_array &rows, const ml_int_array &cols);
        ml_matrix &delete_rows_cols(int row, const ml_int_array &cols);
        ml_matrix &delete_rows_cols(const ml_int_array &rows, int col);
        ml_matrix &delete_rows_cols(int row, int col);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, int n_rows, int *rows);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, const ml_int_array &rows);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, int n_rows, int *rows, int n_cols, int *cols);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, const ml_int_array &rows, const ml_int_array &cols);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, int row, int n_cols, int *cols);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, int row, const ml_int_array &cols);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, int n_rows, int *rows, int col);
        ml_matrix &inc_sub_matrix(const ml_matrix &m, const ml_int_array &rows, int col);
        ml_matrix &inc_rows(const ml_matrix &a, int from_row, int to_row);
        ml_matrix &inc_cols(const ml_matrix &m, int from_col, int to_col);
        ml_matrix &inc_elements(const ml_matrix &values,const ml_matrix &elements);
        ml_matrix &dup_rows(int num);
        ml_matrix &dup_cols(int num);
        ml_matrix &append(const ml_matrix &b);
        ml_matrix &stack(const ml_matrix &b);

        #if 0
        #pragma mark BLAS and LAPACK Utility Functions
        #endif



        double *raw_data() const
        {
            return mat;
        }

        int data_stride() const
        {
            return ml_major_stride(real_rows, real_cols);
        }

        #if 0
        #pragma mark Inspection/Assignment Members
        #endif



        int error;
        int is_empty() const;
        int operator==(const ml_matrix &m) const;
        int operator!=(const ml_matrix &m) const;
        ml_matrix find(ml_comparison op, double target) const;
        void sort(int direction = 0);
        ml_matrix sort_indices(int direction = 0) const;
        void permute(const ml_matrix &elements);
        ml_matrix sum_rows() const;
        ml_matrix sum_cols() const;
        ml_matrix max() const;
        ml_matrix row_max() const;
        ml_matrix col_max() const;
        ml_matrix min() const;
        ml_matrix row_min() const;
        ml_matrix col_min() const;
        ml_matrix interpolate(int col, double value) const;

        inline int rows() const
        {
            return num_rows;
        }

        inline int cols() const
        {
            return num_cols;
        }
        ml_matrix &operator=(const ml_matrix &m);
        ml_matrix &operator=(const char *s);
        void zero();
        void ones();
        void set_val(double v);
        void set_val(const ml_matrix &elements,double v);
        ml_matrix &apply(double (*func)(double));
        void apply_to_elements(const ml_matrix &elements, double (*func)(double));
        void identity();
        void rand_seed(double seed);
        void rand();
        void randn(double std_dev);

        #if 0
        #pragma mark Arithmetic Members
        #endif


        ml_matrix &operator*=(const ml_matrix &m);
        ml_matrix &operator*=(double a);
        ml_matrix &dotstar_eq(const ml_matrix &m);
        ml_matrix &operator/=(const ml_matrix &m);
        ml_matrix &operator/=(double a);
        ml_matrix &dotslash_eq(const ml_matrix &m);
        ml_matrix &operator+=(const ml_matrix &m);
        ml_matrix &operator+=(double a);
        ml_matrix &operator-=(const ml_matrix &m);
        ml_matrix &operator-=(double a);

        #if 0
        #pragma mark Algebraic Members
        #endif


        ml_matrix &abs();
        ml_matrix &ceil();
        ml_matrix &floor();
        ml_matrix &round();
        ml_matrix &rint();
        ml_matrix &dotpow(double n);
        ml_matrix &sqrt();
        ml_matrix &exp();
        ml_matrix &expm1();
        ml_matrix &log();
        ml_matrix &log10();
        ml_matrix &log1p();
        ml_matrix &log2();
        ml_matrix &logb();

        #if 0
        #pragma mark Linear Alegbra Members
        #endif


        ml_matrix &transpose();
        ml_matrix &inv();
        ml_matrix invp(const ml_matrix &n, int k) const;
        ml_matrix &mpow(int n);
        ml_matrix diag() const;
        double trace() const;
        ml_matrix mag() const;
        double vmag() const;
        ml_matrix &cross(const ml_matrix &n);
        double dot(const ml_matrix &n) const;
        ml_matrix skew_symm() const;
        ml_matrix &unit();
        ml_matrix &sign();
        ml_matrix &vunit();
        double norm() const;
        ml_matrix &solve_ax_eq_b(const ml_matrix &a);
        inline ml_matrix &row_mult(int row, double fac)
        {
            double *m2 = ml_matloc(mat, row-1, 0, real_rows, real_cols);
            for (int i = 0; i < num_cols; i++)
            {
                *m2 *= fac;
                ml_incr_col(m2,real_rows, real_cols);
            }
            return *this;
        }

        #if 0
        #pragma mark Trigonometric Members
        #endif


        ml_matrix &sin();
        ml_matrix &cos();
        ml_matrix &tan();
        ml_matrix &csc();
        ml_matrix &sec();
        ml_matrix &cot();
        ml_matrix &asin();
        ml_matrix &acos();
        ml_matrix &atan();
        ml_matrix &acsc();
        ml_matrix &asec();
        ml_matrix &acot();
        ml_matrix &sinh();
        ml_matrix &cosh();
        ml_matrix &tanh();
        ml_matrix &csch();
        ml_matrix &sech();
        ml_matrix &coth();
        ml_matrix &asinh();
        ml_matrix &acosh();
        ml_matrix &atanh();
        ml_matrix &acsch();
        ml_matrix &asech();
        ml_matrix &acoth();
        ml_matrix &atan2(const ml_matrix &m);

    protected:
        // The array containing the elements of the matrix.
        double *mat;
        // The number of rows in the matrix.
        int num_rows;
        // The number of columns in the matrix.
        int num_cols;
        // The number of rows allocated for storage.
        int real_rows;
        // The number of columns allocated for storage.
        int real_cols;
        // The size of the storage array, equal to the product of real_rows and real_cols.
        int len;
        // Character buffer used for to_string() operations.
        char *buff;
        // Size of the character buffer.
        int buffsize;
        // Seed for randomizing operations.
        double seed;
        // Built-in static array for optimizing small matrices.
        double base_mat[ML_ARR_INCR*ML_ARR_INCR];
        // Flag indicating whether the built-in or dynamically allocated array is used.
        int using_base;
        inline void swap_rows(int a, int b);
        inline void row_add(int a, int b);
        inline void row_sub(int k, int j);
};

#if 0
#pragma mark Construction and I/O Non-Members
#endif


ml_matrix linspace(double start, double finish, int n);
ml_matrix colon(double start, double incr, double end);

#if 0
#pragma mark Matrix Manipulation Non-Members
#endif


ml_matrix append_identity(const ml_matrix &b);
ml_matrix dup_rows(const ml_matrix &m, int num);
ml_matrix dup_cols(const ml_matrix &m, int num);
ml_matrix append(const ml_matrix &a, const ml_matrix &b);
ml_matrix stack(const ml_matrix &a, const ml_matrix &b);

#if 0
#pragma mark Comparison, Inspection, Assignment Non-Members
#endif


ml_matrix sort(const ml_matrix &m, int direction = 0);
ml_matrix apply(const ml_matrix &m, double (*func)(double));

#if 0
#pragma mark Arithmetic Non-Members
#endif


ml_matrix operator*(double b,const ml_matrix &a);
ml_matrix operator*(const ml_matrix &a, double b);
ml_matrix operator*(const ml_matrix &a, const ml_matrix &b);
ml_matrix dotstar(const ml_matrix &a, const ml_matrix &b);
ml_matrix operator/(const ml_matrix &a, double b);
ml_matrix operator/(const ml_matrix &a, const ml_matrix &b);
ml_matrix backslash(const ml_matrix &a, const ml_matrix &b);
ml_matrix dotslash(const ml_matrix &a, const ml_matrix &b);
ml_matrix operator+(double b,const ml_matrix &a);
ml_matrix operator+(const ml_matrix &a, double b);
ml_matrix operator+(const ml_matrix &a, const ml_matrix &b);
ml_matrix operator-(double b,const ml_matrix &a);
ml_matrix operator-(const ml_matrix &a, double b);
ml_matrix operator-(const ml_matrix &a, const ml_matrix &b);
ml_matrix operator-(const ml_matrix &a);

#if 0
#pragma mark Algebraic Non-Members
#endif


ml_matrix abs(const ml_matrix &x);
ml_matrix ceil(const ml_matrix &x);
ml_matrix floor(const ml_matrix &x);
ml_matrix round(const ml_matrix &x);
ml_matrix rint(const ml_matrix &x);
ml_matrix dotpow(const ml_matrix&, double n);
ml_matrix sqrt(const ml_matrix &x);
ml_matrix exp(const ml_matrix &x);
ml_matrix expm1(const ml_matrix &x);
ml_matrix log(const ml_matrix &x);
ml_matrix log10(const ml_matrix &x);
ml_matrix log1p(const ml_matrix &x);
ml_matrix log2(const ml_matrix &x);
ml_matrix logb(const ml_matrix &x);

#if 0
#pragma mark Linear Algebra Non-Members
#endif


ml_matrix transpose(const ml_matrix &m);
ml_matrix inv(const ml_matrix &m);
ml_matrix mpow(const ml_matrix &m, int n);
ml_matrix cross(const ml_matrix &m, const ml_matrix &y);
ml_matrix unit(const ml_matrix &m);
ml_matrix sign(const ml_matrix &m);
ml_matrix svd(ml_matrix &u, ml_matrix &v, const ml_matrix &a);
ml_matrix vunit(const ml_matrix &m);
ml_matrix solve_ax_eq_b(const ml_matrix &a,const ml_matrix &b);
bool simplex(ml_matrix &x, ml_matrix c, ml_matrix a, ml_matrix b, double max_val);
void simplex2(ml_matrix a, ml_matrix c, ml_matrix &x, ml_int_array &s, ml_int_array &g, ml_int_array &gu, ml_matrix &b_inv, double max_val);

#if 0
#pragma mark Trigonometric Non-Members
#endif


ml_matrix sin(const ml_matrix &x);
ml_matrix cos(const ml_matrix &x);
ml_matrix tan(const ml_matrix &x);
ml_matrix csc(const ml_matrix &x);
ml_matrix sec(const ml_matrix &x);
ml_matrix cot(const ml_matrix &x);
ml_matrix asin(const ml_matrix &x);
ml_matrix acos(const ml_matrix &x);
ml_matrix atan(const ml_matrix &x);
ml_matrix acsc(const ml_matrix &x);
ml_matrix asec(const ml_matrix &x);
ml_matrix acot(const ml_matrix &x);
ml_matrix atan2(const ml_matrix &x, const ml_matrix &y);
ml_matrix sinh(const ml_matrix &x);
ml_matrix cosh(const ml_matrix &x);
ml_matrix tanh(const ml_matrix &x);
ml_matrix csch(const ml_matrix &x);
ml_matrix sech(const ml_matrix &x);
ml_matrix coth(const ml_matrix &x);
ml_matrix asinh(const ml_matrix &x);
ml_matrix acosh(const ml_matrix &x);
ml_matrix atanh(const ml_matrix &x);
ml_matrix acsch(const ml_matrix &x);
ml_matrix asech(const ml_matrix &x);
ml_matrix acoth(const ml_matrix &x);

}

#endif