#ifndef __MATRIX_H_
#define __MATRIX_H_

#include <vector>

using namespace std;

double **multiply(double **A, int am, int an, double **B, int bm, int bn);
double **inverse(double **A, int N, double &determ);
double **transpose(double **A, int m, int n);

double **mplus(double **A, double **B, int m, int n);
double **mminus(double **A, double **B, int m, int n);

double **sub_matrix(double **A, int m, int n, const vector<int> & row_index, const vector<int> & col_index);

void clear_matrix(double **A, int m, int n);
void display_matrix(double **A, int m, int n);

int cholesky(double **orig, int N, double **aug, int mcol, double **chol, double **cholaug, int ofs);
double **self_colmultiply(double **A, int am, int m_start, int m_end, int n_start, int n_end);
double **quadratic_product(double **A, int am, int an, double **B, int bm, int bn, double **C, int cm, int cn);  // A^T B C

#endif