#include <stdint.h>
#include <math.h>
#include <full_util.h>
template <class real, class integer>
real ** allocate_matrix(integer m, integer n) {
real **mat;
mat=new real *[m];
mat[0]=new real[m*n];
for (integer i=1; i<m; i++) {
mat[i]=mat[0]+n*i;
}
return mat;
}
template <class real, class integer>
void zero_matrix(real ** mat, integer m, integer n) {
for (integer i=0; i<m*n; i++) mat[0][i]=0;
}
template <class real, class integer>
real ** zero_matrix(integer m, integer n) {
real **result;
result=allocate_matrix<real, integer>(m, n);
zero_matrix(result, m, n);
return result;
}
template <class real, class integer>
void identity_matrix(real ** mat, integer m, integer n) {
integer dmin;
zero_matrix(mat, m, n);
if (m < n) dmin=m; else dmin=n;
for (integer i=0; i<dmin; i++) mat[i][i]=1;
}
template <class real, class integer>
real ** identity_matrix(integer m, integer n) {
real **result;
result=allocate_matrix<real, integer>(m, n);
identity_matrix(result, m, n);
return result;
}
template <class real>
void delete_matrix(real ** mat) {
if (mat==NULL) return;
delete[] mat[0];
delete[] mat;
}
template <class real, class integer>
void copy_matrix(real **m1, real **m2, integer m, integer n) {
for (integer i=0; i<m; i++) {
for (integer j=0; j<n; j++) {
m2[i][j]=m1[i][j];
}
}
}
template <class real, class integer>
real ** copy_matrix(real **mat, integer m, integer n) {
real **result;
//this is really fucking annoying, I have to add all these fucking
//template parameters in order for it to compile:
result=allocate_matrix<real, integer>(m, n);
copy_matrix(mat, result, m, n);
return result;
}
//should be more efficient:
template <class real, class integer>
void matrix_mult_t(real **plier, real **cand, real **result, integer m, integer p, integer n) {
for (integer i=0; i<m; i++) {
for (integer j=0; j<n; j++) {
result[i][j]=0;
for (integer k=0; k<p; k++) {
result[i][j]+=plier[k][i]*cand[k][j];
}
}
}
}
template <class real, class integer>
void matrix_mult(real **plier, real **cand, real **result, integer m, integer p, integer n) {
for (integer i=0; i<m; i++) {
for (integer j=0; j<n; j++) {
result[i][j]=0;
for (integer k=0; k<p; k++) {
result[i][j]+=plier[i][k]*cand[k][j];
}
}
}
}
template <class real, class integer>
real ** matrix_mult(real **plier, real **cand, integer m, integer p, integer n) {
real **result;
result=allocate_matrix<real, integer>(m, n);
matrix_mult(plier, cand, result, m, p, n);
return result;
}
template <class real, class integer>
void vector_mult(real **plier, real *cand, real *result, integer m, integer n) {
for (integer i=0; i<m; i++) {
result[i]=0;
for (integer j=0; j<n; j++) result[i]+=plier[i][j]*cand[j];
}
}
template <class real, class integer>
real * vector_mult(real **plier, real *cand, integer m, integer n) {
real *result;
result=new real[n];
vector_mult(plier, cand, result, m, n);
return result;
}
template <class real, class integer>
void left_vec_mult(real *plier, real **cand, real *result, integer m, integer n) {
for (integer i=0; i<n; i++) {
result[i]=0;
for (integer j=0; j<m; j++) result[i]+=plier[j]*cand[j][i];
}
}
template <class real, class integer>
real * left_vec_mult(real *plier, real **cand, integer m, integer n) {
real *result;
result=new real[n];
left_vec_mult<real, integer>(plier, cand, result, m, n);
return result;
}
template <class real, class integer>
void matrix_add(real **mat1, real ** mat2, integer m, integer n) {
for (integer i=0; i<m; i++) {
for (integer j=0; j<n; j++) {
mat1[i][j]=mat1[i][j]+mat2[i][j];
}
}
}
//for square matrices (inplace):
template <class real, class integer>
void matrix_transpose(real **mat, integer m) {
real swap;
for (integer i=0; i<m; i++) {
for (integer j=0; j<i; j++) {
swap=mat[j][i];
mat[j][i]=mat[i][j];
mat[i][j]=swap;
}
}
}
#define MAXLL 5000
template <class real, class integer>
real ** matrix_transpose(real **mat, integer m, integer n) {
real **matnew;
//fuuuuuuuuuuuu
matnew=allocate_matrix<real, integer>(n, m);
for (integer i=0; i<m; i++) {
for (integer j=0; j<n; j++) {
matnew[j][i]=mat[i][j];
}
}
return matnew;
}
#define MAXLL 5000
template <class real, class integer>
real ** scan_matrix(FILE *fptr, integer &m, integer &n) {
char line[MAXLL];
real **mat;
int line_pos=0;
int lc;
if (fgets(line, MAXLL, fptr)==NULL) return NULL;
sscanf(line, "%d %d\n", &m, &n);
//fuuuuuu fufufuuuuuuuuuu..... cccccck
mat=allocate_matrix<real, integer>(m, n);
for (integer i=0; i<m; i++) {
if (fgets(line, MAXLL, fptr)==NULL) {
fprintf(stderr, "scan_matrix: file ended before finished scanning\n");
break;
}
for (integer j=0; j<n; j++) {
sscanf(line+line_pos, "%f%n", &mat[i][j], &lc);
line_pos+=lc;
}
}
return mat;
}
template <class real, class integer>
void print_matrix(FILE *fptr, real **mat, integer m, integer n) {
fprintf(fptr, "%d %d\n", m, n);
for (integer i=0; i<m; i++) {
for (integer j=0; j<n; j++) {
fprintf(fptr, "%f ", mat[i][j]);
}
fprintf(fptr, "\n");
}
}
//getting too long to be inlined...
template <class real, class integer>
real ** read_matrix(FILE *fptr, integer &m, integer &n) {
real **mat;
long check;
fread(&n, 1, sizeof(n), fptr);
check=ftell(fptr);
//not a seekable stream--use a linked list to get the data:
if (check < 0) {
linked_list<real_a> data;
real_a *data2;
real_a val=0;
long ntot=0;
while (fread(&val, sizeof(val), 1, fs)==1) data.add(val);
data2=data.make_array(ntot);
//assert(ntot % nvar==0);
ntrain=ntot/nvar;
train=new real_a *[ntrain];
train[0]=data2;
for (nel_ta i=0; i<ntrain; i++) train[i]=train[0]+i*nvar;
if (ntot % n != 0) m=-1;
} else {
fseek(fptr, 0, SEEK_END);
m=(ftell(fptr)-sizeof(n))/sizeof(real)/n;
fseek(fptr, sizeof(n), SEEK_SET);
//printf("read_matrix: reading a [%dx%d] matrix\n", m, n);
mat=allocate_matrix<real, integer>(m, n);
fread(mat[0], sizeof(real), m*n, fptr);
}
return mat;
}
template <class real, class integer>
size_t write_matrix(FILE *fptr, real **mat, integer m, integer n) {
size_t nwrit;
nwrit+=fwrite(&n, 1, sizeof(n), fptr);
nwrit+=fwrite(mat[0], sizeof(real), m*n, fptr);
return nwrit;
}
template <class real, class integer>
real matrix_norm(real **mat, integer m, integer n) {
real result=0;
for (integer i=0; i<m*n; i++) result+=mat[0][i]*mat[0][i];
return sqrt(result);
}
template float ** zero_matrix<float, int16_t>(int16_t, int16_t);
template float ** zero_matrix<float, int32_t>(int32_t, int32_t);
template float ** zero_matrix<float, int64_t>(int16_t, int64_t);
template double ** zero_matrix<double, int16_t>(int16_t, int16_t);
template double ** zero_matrix<double, int32_t>(int32_t, int32_t);
template double ** zero_matrix<double, int64_t>(int16_t, int64_t);
template float ** identity_matrix<float, int16_t>(int16_t, int16_t);
template float ** identity_matrix<float, int32_t>(int32_t, int32_t);
template float ** identity_matrix<float, int64_t>(int16_t, int64_t);
template double ** identity_matrix<double, int16_t>(int16_t, int16_t);
template double ** identity_matrix<double, int32_t>(int32_t, int32_t);
template double ** identity_matrix<double, int64_t>(int16_t, int64_t);
template void delete_matrix<float>(float **);
template void delete_matrix<double>(double **);
template float ** copy_matrix<float, int16_t>(float **, int16_t, int16_t);
template float ** copy_matrix<float, int32_t>(float **, int32_t, int32_t);
template float ** copy_matrix<float, int64_t>(float **, int16_t, int64_t);
template double ** copy_matrix<double, int16_t>(double **, int16_t, int16_t);
template double ** copy_matrix<double, int32_t>(double **, int32_t, int32_t);
template double ** copy_matrix<double, int64_t>(double **, int64_t);
template float ** matrix_mult<float, int16_t>(float **, float **, int16_t, int16_t, int16_t);
template float ** matrix_mult<float, int32_t>(float **, float **, int32_t, int32_t, int32_t);
template float ** matrix_mult<float, int64_t>(float **, float **, int64_t, int64_t, int64_t);
template double ** matrix_mult<double, int16_t>(double **, double **, int16_t, int16_t, int16_t);
template double ** matrix_mult<double, int32_t>(double **, double **, int32_t, int32_t, int32_t);
template double ** matrix_mult<double, int64_t>(double **, double **, int64_t, int64_t, int64_t);
template void matrix_mult_t<float, int16_t>(float **, float **, float **, int16_t, int16_t, int16_t);
template void matrix_mult_t<float, int32_t>(float **, float **, float **, int32_t, int32_t, int32_t);
template void matrix_mult_t<float, int64_t>(float **, float **, float **, int64_t, int64_t, int64_t);
template void matrix_mult_t<double, int16_t>(double **, double **, double **,
int16_t, int16_t, int16_t);
template void matrix_mult_t<double, int32_t>(double **, double **, double **,
int32_t, int32_t, int32_t);
template void matrix_mult_t<double, int64_t>(double **, double **, double **,
int64_t, int64_t, int64_t);
template float * vector_mult<float, int16_t>(float **, float *, int16_t, int16_t);
template float * vector_mult<float, int32_t>(float **, float *, int32_t, int32_t);
template float * vector_mult<float, int64_t>(float **, float *, int64_t, int64_t);
template double * vector_mult<double, int16_t>(double **, double *, int16_t, int16_t);
template double * vector_mult<double, int32_t>(double **, double *, int32_t, int32_t);
template double * vector_mult<double, int64_t>(double **, double *, int64_t, int64_t);
template float * left_vec_mult<float, int16_t>(float *, float **, int16_t, int16_t);
template float * left_vec_mult<float, int32_t>(float *, float **, int32_t, int32_t);
template float * left_vec_mult<float, int64_t>(float *, float **, int64_t, int64_t);
template double * left_vec_mult<double, int16_t>(double *, double **, int16_t, int16_t);
template double * left_vec_mult<double, int32_t>(double *, double **, int32_t, int32_t);
template double * left_vec_mult<double, int64_t>(double *, double **, int64_t, int64_t);
//fuck, half of this file is just gonna be template instantiations...
template float ** matrix_transpose<float, int16_t>(float **, int16_t, int16_t);
template float ** matrix_transpose<float, int32_t>(float **, int32_t, int32_t);
template float ** matrix_transpose<float, int64_t>(float **, int64_t, int64_t);
template double ** matrix_transpose<double, int16_t>(double **, int16_t, int16_t);
template double ** matrix_transpose<double, int32_t>(double **, int32_t, int32_t);
template double ** matrix_transpose<double, int64_t>(double **, int64_t, int64_t);
template float ** scan_matrix<float, int16_t>(FILE *fs, int16_t &, int16_t &);
template float ** scan_matrix<float, int32_t>(FILE *fs, int32_t &, int32_t &);
template float ** scan_matrix<float, int64_t>(FILE *fs, int64_t &, int64_t &);
template double ** scan_matrix<double, int16_t>(FILE *fs, int16_t, int16_t);
template double ** scan_matrix<double, int32_t>(FILE *fs, int32_t, int32_t);
template double ** scan_matrix<double, int64_t>(FILE *fs, int64_t, int64_t);
template float ** read_matrix<float, int16_t>(FILE *fs, int16_t &, int16_t &);
template float ** read_matrix<float, int32_t>(FILE *fs, int32_t &, int32_t &);
template float ** read_matrix<float, int64_t>(FILE *fs, int64_t &, int64_t &);
template double ** read_matrix<double, int16_t>(FILE *fs, int16_t, int16_t);
template double ** read_matrix<double, int32_t>(FILE *fs, int32_t, int32_t);
template double ** read_matrix<double, int64_t>(FILE *fs, int64_t, int64_t);
template void print_matrix<float, int16_t>(FILE *fs, float **, int16_t &, int16_t &);
template void print_matrix<float, int32_t>(FILE *fs, float **, int32_t &, int32_t &);
template void print_matrix<float, int64_t>(FILE *fs, int64_t &, int64_t &);
template void print_matrix<double, int16_t>(FILE *fs, double **, int16_t, int16_t);
template void print_matrix<double, int32_t>(FILE *fs, double **, int32_t, int32_t);
template void print_matrix<double, int64_t>(FILE *fs, double **, int64_t, int64_t);
template size_t write_matrix<float, int16_t>(FILE *fs, float **, int16_t &, int16_t &);
template size_t write_matrix<float, int32_t>(FILE *fs, float **, int32_t &, int32_t &);
template size_t write_matrix<float, int64_t>(FILE *fs, int64_t &, int64_t &);
template size_t write_matrix<double, int16_t>(FILE *fs, double **, int16_t, int16_t);
template size_t write_matrix<double, int32_t>(FILE *fs, double **, int32_t, int32_t);
template size_t write_matrix<double, int64_t>(FILE *fs, double **, int64_t, int64_t);
template float matrix_norm<float, int16_t>(float **, int16_t, int16_t);
template float matrix_norm<float, int32_t>(float **, int32_t, int32_t);
template float matrix_norm<float, int64_t>(float **, int64_t, int64_t);
template double matrix_norm<double, int16_t>(double **, int16_t, int16_t);
template double matrix_norm<double, int32_t>(double **, int32_t, int32_t);
template double matrix_norm<double, int64_t>(double **, int64_t, int64_t);