https://github.com/stamatak/standard-RAxML
Tip revision: a3d300633918f16910238f90924d1cfe09c59b13 authored by stamatak on 23 March 2016, 09:41:50 UTC
integrated Sarah's code for randomly subsampling quartets without replacement.
integrated Sarah's code for randomly subsampling quartets without replacement.
Tip revision: a3d3006
eigen.c
/* RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
* Copyright August 2006 by Alexandros Stamatakis
*
* Partially derived from
* fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
*
* and
*
* Programs of the PHYLIP package by Joe Felsenstein.
*
* This program is free software; you may redistribute it and/or modify its
* under the terms of the GNU General Public License as published by the Free
* Software Foundation; either version 2 of the License, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
* for more details.
*
*
* For any other enquiries send an Email to Alexandros Stamatakis
* Alexandros.Stamatakis@epfl.ch
*
* When publishing work that is based on the results from RAxML-VI-HPC please cite:
*
* Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands
* of taxa and mixed models".
* Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
*/
#include <math.h>
#include <assert.h>
static void mytred2(double **a, const int n, double *d, double *e)
{
int l, k, j, i;
double scale, hh, h, g, f;
for (i = n; i > 1; i--)
{
l = i - 1;
h = 0.0;
scale = 0.0;
if (l > 1)
{
for (k = 1; k <= l; k++)
scale += fabs(a[k - 1][i - 1]);
if (scale == 0.0)
e[i - 1] = a[l - 1][i - 1];
else
{
for (k = 1; k <= l; k++)
{
a[k - 1][i - 1] /= scale;
h += a[k - 1][i - 1] * a[k - 1][i - 1];
}
f = a[l - 1][i - 1];
g = ((f > 0) ? -sqrt(h) : sqrt(h)); /* diff */
e[i - 1] = scale * g;
h -= f * g;
a[l - 1][i - 1] = f - g;
f = 0.0;
for (j = 1; j <= l; j++)
{
a[i - 1][j - 1] = a[j - 1][i - 1] / h;
g = 0.0;
for (k = 1; k <= j; k++)
g += a[k - 1][j - 1] * a[k - 1][i - 1];
for (k = j + 1; k <= l; k++)
g += a[j - 1][k - 1] * a[k - 1][i - 1];
e[j - 1] = g / h;
f += e[j - 1] * a[j - 1][i - 1];
}
hh = f / (h + h);
for (j = 1; j <= l; j++)
{
f = a[j - 1][i - 1];
g = e[j - 1] - hh * f;
e[j - 1] = g;
for (k = 1; k <= j; k++)
a[k - 1][j - 1] -= (f * e[k - 1] + g * a[k - 1][i - 1]);
}
}
}
else
e[i - 1] = a[l - 1][i - 1];
d[i - 1] = h;
}
d[0] = 0.0;
e[0] = 0.0;
for (i = 1; i <= n; i++)
{
l = i - 1;
if (d[i - 1] != 0.0)
{
for (j = 1; j <= l; j++)
{
g = 0.0;
for (k = 1; k <= l; k++)
g += a[k - 1][i - 1] * a[j - 1][k - 1];
for(k = 1; k <= l; k++)
a[j - 1][k - 1] -= g * a[i - 1][k - 1];
}
}
d[i - 1] = a[i - 1][i - 1];
a[i - 1][i - 1] = 1.0;
for (j = 1; j <= l; j++)
a[i - 1][j - 1] = a[j - 1][i - 1] = 0.0;
}
}
/*#define MYSIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))*/
static int mytqli(double *d, double *e, const int n, double **z)
{
int m, l, iter, i, k;
double s, r, p, g, f, dd, c, b;
for (i = 2; i <= n; i++)
e[i - 2] = e[i - 1];
e[n - 1] = 0.0;
for (l = 1; l <= n; l++)
{
iter = 0;
do
{
for (m = l; m <= n - 1; m++)
{
dd = fabs(d[m - 1]) + fabs(d[m]);
if (fabs(e[m - 1]) + dd == dd)
break;
}
if (m != l)
{
assert(iter < 30);
g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
r = sqrt((g * g) + 1.0);
g = d[m - 1] - d[l - 1] + e[l - 1] / (g + ((g < 0)?-fabs(r):fabs(r)));/*MYSIGN(r, g));*/
s = c = 1.0;
p = 0.0;
for (i = m - 1; i >= l; i--)
{
f = s * e[i - 1];
b = c * e[i - 1];
if (fabs(f) >= fabs(g))
{
c = g / f;
r = sqrt((c * c) + 1.0);
e[i] = f * r;
c *= (s = 1.0 / r);
}
else
{
s = f / g;
r = sqrt((s * s) + 1.0);
e[i] = g * r;
s *= (c = 1.0 / r);
}
g = d[i] - p;
r = (d[i - 1] - g) * s + 2.0 * c * b;
p = s * r;
d[i] = g + p;
g = c * r - b;
for (k = 1; k <= n; k++)
{
f = z[i][k-1];
z[i][k-1] = s * z[i - 1][k - 1] + c * f;
z[i - 1][k - 1] = c * z[i - 1][k - 1] - s * f;
}
}
d[l - 1] = d[l - 1] - p;
e[l - 1] = g;
e[m - 1] = 0.0;
}
}
while (m != l);
}
return (1);
}
void makeEigen(double **_a, const int n, double *d, double *e)
{
mytred2(_a, n, d, e);
mytqli(d, e, n, _a);
}