/* MIPoD 1.0 basic matrix operations header file */ /* copyright 2007 Paul Hohenlohe */ /* permission is granted to use and modify this software _provided that_ any publication resulting from use of this software or any part or modification thereof contains at least one of the following citations: Hohenlohe, P.A. 2007. MIPoD: Microevolutionary Inference from Patterns of Divergence. Software available at http://oregonstate.edu/~hohenlop/software Hohenlohe, P.A. & S.J. Arnold. 2008. MIPoD: A Hypothesis-Testing Framework for Microevolutionary Inference from Patterns of Divergence. American Naturalist 171(3): 366-385. questions/comments: hohenlop@science.oregonstate.edu */ /* Functions defined herein: minimum(x,y) maximum(x,y) kroneckerdelta(x,y) bn logdeterminant(double **a, int n) int SquareKroneckerProduct(double **a, double **b, double **c, int n, int m) int LUinverse(double **a, double **b, int n) int linearsolve(double **mat, double *x, double *y, int n) int rowechelon(double **mat, int r, int c) double chisquare(double x, int d) int eigen(double **mat, double **out, int n) int inveigen(double **evecs, double **out, int n, int y) int ordereigen(double **evecs, int n) */ #define minimum(x,y) ( ((x) < (y)) ? (x) : (y) ) #define maximum(x,y) ( ((x) > (y)) ? (x) : (y) ) #define kroneckerdelta(x,y) ( ((x) == (y)) ? (1) : (0) ) typedef struct { int i; double d; } bn; /* store big numbers as i*e^d, so d = ln(|bn|) and i is the sign */ bn logdeterminant(double **a, int n) { /* uses PLU decomposition algorithm from Demmel 1997; returns determinant of a as a bn */ int i, j, k, sign, zero, row, col; double m, tmp, **b; bn det; if ((b = malloc((n+1)*sizeof(double *))) == NULL) { printf ("memory error"); exit(0); } for (i=0; i<(n+1); i++) if ((b[i] = malloc((n+1)*sizeof(double))) == NULL) { printf ("memory error"); exit(0); } for(i=0; i m) { m = tmp; row = k; col = j; } } } if(m == 0) { zero = 1; continue; } if(row != i) { for(j=0; j 0) det.i = 1; else det.i = -1; } for (i=0; i<(n+1); i++) free(b[i]); free (b); return (det); } int SquareKroneckerProduct(double **a, double **b, double **c, int n, int m) { /* Populates square matrix c (size (n*m)x(n*m)) with the Kronecker product of square matrices a (size nxn) and b (size mxm) */ int i,j,k,l; for(i=0; i m) {m = tmp; row = k; col = j; } } } if(m == 0) {zero = 1; continue; } if(row != i) { /*switch rows */ for(k=0; k=0; i--) { /* invert upper factor */ UL[i][i] = 1/b[i][i]; for(j=(i+1); j 1) { /* If df is odd and >1, sum forms the kernel of larger exact calc. */ for(q=0.5; q fabs(tmp)) { k = i; l = j; tmp = tmpmat[i][j]; } } } while (fabs(tmp) > 0.00001*sqrt(fabs(tmpmat[k][k]*tmpmat[l][l]))) { if (tmpmat[k][k] == tmpmat[l][l]) rot = 0.25*M_PI; else rot = 0.5*atan(2*tmpmat[k][l]/(tmpmat[k][k]-tmpmat[l][l])); c = cos(rot); s = sin(rot); /* lower triangle is rotated, then transfered to upper; k fabs(tmp)) { k = i; l = j; tmp = tmpmat[i][j]; } } } } for(i=0; i tmp) { tmp = fabs(tmpmat[j][j]); k = j; } } out[n][i] = tmpmat[k][k]; tmpmat[k][k] = 0; } for(i=0; i tmp) { tmp = fabs(tmpevecs[n][j]); k = j; } } for(j=0; j<(n+1); j++) evecs[j][i] = tmpevecs[j][k]; tmpevecs[n][k] = 0; } for(i=0; i<(n+1); i++) free(tmpevecs[i]); free(tmpevecs); return 0; }