source: tags/initial/GDE/MOLPHY/tranprb.c

Last change on this file was 2, checked in by oldcode, 24 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.3 KB
Line 
1/*
2 * tranprob.c   Adachi, J.   1995.10.19
3 * Copyright (C) 1992-1995 J. Adachi & M. Hasegawa; All rights reserved.
4 */
5
6#define SH 1 /* suggested by K.Strimmer & A.v.Haeseler */
7
8#define TPMWRITE 0
9
10#include "protml.h"
11
12
13void 
14elmhes(a, ordr, n)
15        dmattpmty a;
16        int *ordr;
17        int n;
18{
19        int m, j, i;
20        double y, x;
21
22        for (i = 0; i < n; i++)
23                ordr[i] = 0;
24        for (m = 2; m < n; m++) {
25                x = 0.0;
26                i = m;
27                for (j = m; j <= n; j++) {
28                        if (fabs(a[j - 1][m - 2]) > fabs(x)) {
29                                x = a[j - 1][m - 2];
30                                i = j;
31                        }
32                }
33                ordr[m - 1] = i;      /* vector */
34                if (i != m) {
35                        for (j = m - 2; j < n; j++) {
36                                y = a[i - 1][j];
37                                a[i - 1][j] = a[m - 1][j];
38                                a[m - 1][j] = y;
39                        }
40                        for (j = 0; j < n; j++) {
41                                y = a[j][i - 1];
42                                a[j][i - 1] = a[j][m - 1];
43                                a[j][m - 1] = y;
44                        }
45                }
46                if (x != 0.0) {
47                        for (i = m; i < n; i++) {
48                                y = a[i][m - 2];
49                                if (y != 0.0) {
50                                        y /= x;
51                                        a[i][m - 2] = y;
52                                        for (j = m - 1; j < n; j++)
53                                                a[i][j] -= y * a[m - 1][j];
54                                        for (j = 0; j < n; j++)
55                                                a[j][m - 1] += y * a[j][i];
56                                }
57                        }
58                }
59        }
60} /*_ elmhes */
61
62
63void 
64eltran(a, zz, ordr, n)
65        dmattpmty a, zz;
66        int *ordr;
67        int n;
68{
69        int i, j, m;
70
71        for (i = 0; i < n; i++) {
72                for (j = i + 1; j < n; j++) {
73                        zz[i][j] = 0.0;
74                        zz[j][i] = 0.0;
75                }
76                zz[i][i] = 1.0;
77        }
78        if (n <= 2)
79                return;
80        for (m = n - 1; m >= 2; m--) {
81                for (i = m; i < n; i++)
82                        zz[i][m - 1] = a[i][m - 2];
83                i = ordr[m - 1];
84                if (i != m) {
85                        for (j = m - 1; j < n; j++) {
86                                zz[m - 1][j] = zz[i - 1][j];
87                                zz[i - 1][j] = 0.0;
88                        }
89                        zz[i - 1][m - 1] = 1.0;
90                }
91        }
92} /*_ eltran */
93
94
95static void 
96cdiv(ar, ai, br, bi, cr, ci)
97        double ar, ai, br, bi, *cr, *ci;
98{
99        double s, ars, ais, brs, bis;
100
101        s = fabs(br) + fabs(bi);
102        ars = ar / s;
103        ais = ai / s;
104        brs = br / s;
105        bis = bi / s;
106        s = brs * brs + bis * bis;
107        *cr = (ars * brs + ais * bis) / s;
108        *ci = (ais * brs - ars * bis) / s;
109} /*_ cdiv */
110
111
112void 
113hqr2(n, low, hgh, err, h, zz, wr, wi)
114        int n, low, hgh, *err;
115        dmattpmty h, zz;
116        double *wr, *wi;
117{
118        int i, j, k, l, m, en, na, itn, its;
119        double p, q, r, s, t, w, x, y, ra, sa, vi, vr, z, norm, tst1, tst2;
120        boolean notlas;
121
122        *err = 0;
123        norm = 0.0;
124        k = 1;
125        /* store roots isolated by balanc and compute matrix norm */
126        for (i = 0; i < n; i++) {
127                for (j = k - 1; j < n; j++)
128                        norm += fabs(h[i][j]);
129                k = i + 1;
130                if (i + 1 < low || i + 1 > hgh) {
131                        wr[i] = h[i][i];
132                        wi[i] = 0.0;
133                }
134        }
135        en = hgh;
136        t = 0.0;
137        itn = n * 30;
138
139        while (en >= low) {            /* search for next eigenvalues */
140                its = 0;
141                na = en - 1;
142
143                while (en >= 1) {      /* infinietr loop */
144                        /* look for single small sub-diagonal element */
145                        for (l = en; l > low; l--) {
146                                s = fabs(h[l - 2][l - 2]) + fabs(h[l - 1][l - 1]);
147                                if (s == 0.0)
148                                        s = norm;
149                                tst1 = s;
150                                tst2 = tst1 + fabs(h[l - 1][l - 2]);
151                                if (tst2 == tst1)
152                                        goto L100;
153                        }
154                        l = low;
155        L100:
156                        x = h[en - 1][en - 1];  /* form shift */
157                        if (l == en || l == na)
158                                break;
159                        if (itn == 0) { /* all eigenvalues have not converged */
160                                *err = en;
161                                goto Lerror;
162                        }
163                        y = h[na - 1][na - 1];
164                        w = h[en - 1][na - 1] * h[na - 1][en - 1];
165                        /* form exceptional shift */
166                        if (its == 10 || its == 20) {
167                                t += x;
168                                for (i = low - 1; i < en; i++)
169                                        h[i][i] -= x;
170                                s = fabs(h[en - 1][na - 1]) + fabs(h[na - 1][en - 3]);
171                                x = 0.75 * s;
172                                y = x;
173                                w = -0.4375 * s * s;
174                        }
175                        its++;
176                        itn--;
177                        /* look for two consecutive small sub-diagonal elements */
178                        for (m = en - 2; m >= l; m--) {
179                                z = h[m - 1][m - 1];
180                                r = x - z;
181                                s = y - z;
182                                p = (r * s - w) / h[m][m - 1] + h[m - 1][m];
183                                q = h[m][m] - z - r - s;
184                                r = h[m + 1][m];
185                                s = fabs(p) + fabs(q) + fabs(r);
186                                p /= s;
187                                q /= s;
188                                r /= s;
189                                if (m == l)
190                                        break;
191                                tst1 = fabs(p) * (fabs(h[m-2][m-2]) + fabs(z) + fabs(h[m][m]));
192                                tst2 = tst1 + fabs(h[m - 1][m - 2]) * (fabs(q) + fabs(r));
193                                if (tst2 == tst1)
194                                        break;
195                        }
196
197                        for (i = m + 2; i <= en; i++) {
198                                h[i - 1][i - 3] = 0.0;
199                                if (i != m + 2)
200                                        h[i - 1][i - 4] = 0.0;
201                        }
202
203                        /* double qr step involving rows l to en and columns m to en */
204                        for (k = m; k <= na; k++) {
205                                notlas = (k != na);
206                                if (k != m) {
207                                        p = h[k - 1][k - 2];
208                                        q = h[k][k - 2];
209                                        r = 0.0;
210                                        if (notlas)
211                                                r = h[k + 1][k - 2];
212                                        x = fabs(p) + fabs(q) + fabs(r);
213                                        if (x != 0.0) {
214                                                p /= x;
215                                                q /= x;
216                                                r /= x;
217                                        }
218                                }
219                                if (x != 0.0) {
220                                        if (p < 0.0) /* sign */
221                                                s = - sqrt(p * p + q * q + r * r);
222                                        else
223                                                s = sqrt(p * p + q * q + r * r);
224                                        if (k != m)
225                                                h[k - 1][k - 2] = -s * x;
226                                        else {
227                                                if (l != m)
228                                                        h[k - 1][k - 2] = -h[k - 1][k - 2];
229                                        }
230                                        p += s;
231                                        x = p / s;
232                                        y = q / s;
233                                        z = r / s;
234                                        q /= p;
235                                        r /= p;
236                                        if (!notlas) {
237                                                for (j = k - 1; j < n; j++) {   /* row modification */
238                                                        p = h[k - 1][j] + q * h[k][j];
239                                                        h[k - 1][j] -= p * x;
240                                                        h[k][j] -= p * y;
241                                                }
242                                                j = (en < (k + 3)) ? en : (k + 3); /* min */
243                                                for (i = 0; i < j; i++) {       /* column modification */
244                                                        p = x * h[i][k - 1] + y * h[i][k];
245                                                        h[i][k - 1] -= p;
246                                                        h[i][k] -= p * q;
247                                                }
248                                                /* accumulate transformations */
249                                                for (i = low - 1; i < hgh; i++) {
250                                                        p = x * zz[i][k - 1] + y * zz[i][k];
251                                                        zz[i][k - 1] -= p;
252                                                        zz[i][k] -= p * q;
253                                                }
254                                        } else {
255                                                for (j = k - 1; j < n; j++) {   /* row modification */
256                                                        p = h[k - 1][j] + q * h[k][j] + r * h[k + 1][j];
257                                                        h[k - 1][j] -= p * x;
258                                                        h[k][j] -= p * y;
259                                                        h[k + 1][j] -= p * z;
260                                                }
261                                                j = (en < (k + 3)) ? en : (k + 3); /* min */
262                                                for (i = 0; i < j; i++) {       /* column modification */
263                                                        p = x * h[i][k - 1] + y * h[i][k] + z * h[i][k + 1];
264                                                        h[i][k - 1] -= p;
265                                                        h[i][k] -= p * q;
266                                                        h[i][k + 1] -= p * r;
267                                                }
268                                                /* accumulate transformations */
269                                                for (i = low - 1; i < hgh; i++) {
270                                                        p = x * zz[i][k-1] + y * zz[i][k] + z * zz[i][k+1];
271                                                        zz[i][k - 1] -= p;
272                                                        zz[i][k] -= p * q;
273                                                        zz[i][k + 1] -= p * r;
274                                                }
275                                        }
276                                }
277                        }              /* for k */
278                }                      /* while infinite loop */
279
280
281                if (l == en) {         /* one root found */
282                        h[en - 1][en - 1] = x + t;
283                        wr[en - 1] = h[en - 1][en - 1];
284                        wi[en - 1] = 0.0;
285                        en = na;
286                        continue;
287                }
288                y = h[na - 1][na - 1];
289                w = h[en - 1][na - 1] * h[na - 1][en - 1];
290                p = (y - x) / 2.0;
291                q = p * p + w;
292                z = sqrt(fabs(q));
293                h[en - 1][en - 1] = x + t;
294                x = h[en - 1][en - 1];
295                h[na - 1][na - 1] = y + t;
296                if (q >= 0.0) {        /* real pair */
297                        if (p < 0.0) /* sign */
298                                z = p - fabs(z);
299                        else
300                                z = p + fabs(z);
301                        wr[na - 1] = x + z;
302                        wr[en - 1] = wr[na - 1];
303                        if (z != 0.0)
304                                wr[en - 1] = x - w / z;
305                        wi[na - 1] = 0.0;
306                        wi[en - 1] = 0.0;
307                        x = h[en - 1][na - 1];
308                        s = fabs(x) + fabs(z);
309                        p = x / s;
310                        q = z / s;
311                        r = sqrt(p * p + q * q);
312                        p /= r;
313                        q /= r;
314                        for (j = na - 1; j < n; j++) {  /* row modification */
315                                z = h[na - 1][j];
316                                h[na - 1][j] = q * z + p * h[en - 1][j];
317                                h[en - 1][j] = q * h[en - 1][j] - p * z;
318                        }
319                        for (i = 0; i < en; i++) {      /* column modification */
320                                z = h[i][na - 1];
321                                h[i][na - 1] = q * z + p * h[i][en - 1];
322                                h[i][en - 1] = q * h[i][en - 1] - p * z;
323                        }
324                        /* accumulate transformations */
325                        for (i = low - 1; i < hgh; i++) {
326                                z = zz[i][na - 1];
327                                zz[i][na - 1] = q * z + p * zz[i][en - 1];
328                                zz[i][en - 1] = q * zz[i][en - 1] - p * z;
329                        }
330                } else {               /* complex pair */
331                        wr[na - 1] = x + p;
332                        wr[en - 1] = x + p;
333                        wi[na - 1] = z;
334                        wi[en - 1] = -z;
335                }
336                en -= 2;
337        }                              /* while en >= low */
338
339        /* backsubstitute to find vectors of upper triangular form */
340        if (norm != 0.0) {
341                for (en = n; en >= 1; en--) {
342                        p = wr[en - 1];
343                        q = wi[en - 1];
344                        na = en - 1;
345                        if (q == 0.0) {/* real vector */
346                                m = en;
347                                h[en - 1][en - 1] = 1.0;
348                                if (na != 0) {
349                                        for (i = en - 2; i >= 0; i--) {
350                                                w = h[i][i] - p;
351                                                r = 0.0;
352                                                for (j = m - 1; j < en; j++)
353                                                        r += h[i][j] * h[j][en - 1];
354                                                if (wi[i] < 0.0) {
355                                                        z = w;
356                                                        s = r;
357                                                } else {
358                                                        m = i + 1;
359                                                        if (wi[i] == 0.0) {
360                                                                t = w;
361                                                                if (t == 0.0) {
362                                                                        tst1 = norm;
363                                                                        t = tst1;
364                                                                        do {
365                                                                                t = 0.01 * t;
366                                                                                tst2 = norm + t;
367                                                                        } while (tst2 > tst1);
368                                                                }
369                                                                h[i][en - 1] = -(r / t);
370                                                        } else {        /* solve real equations */
371                                                                x = h[i][i + 1];
372                                                                y = h[i + 1][i];
373                                                                q = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i];
374                                                                t = (x * s - z * r) / q;
375                                                                h[i][en - 1] = t;
376                                                                if (fabs(x) > fabs(z))
377                                                                        h[i + 1][en - 1] = (-r - w * t) / x;
378                                                                else
379                                                                        h[i + 1][en - 1] = (-s - y * t) / z;
380                                                        }
381                                                        /* overflow control */
382                                                        t = fabs(h[i][en - 1]);
383                                                        if (t != 0.0) {
384                                                                tst1 = t;
385                                                                tst2 = tst1 + 1.0 / tst1;
386                                                                if (tst2 <= tst1) {
387                                                                        for (j = i; j < en; j++)
388                                                                                h[j][en - 1] /= t;
389                                                                }
390                                                        }
391                                                }
392                                        }
393                                }
394                        } else if (q > 0.0) {
395                                m = na;
396                                if (fabs(h[en - 1][na - 1]) > fabs(h[na - 1][en - 1])) {
397                                        h[na - 1][na - 1] = q / h[en - 1][na - 1];
398                                        h[na - 1][en - 1] = (p - h[en-1][en-1]) / h[en-1][na-1];
399                                } else
400                                        cdiv(0.0, -h[na-1][en-1], h[na-1][na-1] - p, q, &h[na-1]
401                                             [na-1], &h[na-1][en-1]);
402                                h[en - 1][na - 1] = 0.0;
403                                h[en - 1][en - 1] = 1.0;
404                                if (en != 2) {
405                                        for (i = en - 3; i >= 0; i--) {
406                                                w = h[i][i] - p;
407                                                ra = 0.0;
408                                                sa = 0.0;
409                                                for (j = m - 1; j < en; j++) {
410                                                        ra += h[i][j] * h[j][na - 1];
411                                                        sa += h[i][j] * h[j][en - 1];
412                                                }
413                                                if (wi[i] < 0.0) {
414                                                        z = w;
415                                                        r = ra;
416                                                        s = sa;
417                                                } else {
418                                                        m = i + 1;
419                                                        if (wi[i] == 0.0)
420                                                                cdiv(-ra, -sa, w, q, &h[i][na-1], &h[i][en-1]);
421                                                        else {  /* solve complex equations */
422                                                                x = h[i][i + 1];
423                                                                y = h[i + 1][i];
424                                                                vr = (wr[i]-p) * (wr[i]-p) + wi[i]*wi[i] - q*q;
425                                                                vi = (wr[i] - p) * 2.0 * q;
426                                                                if (vr == 0.0 && vi == 0.0) {
427                                                                        tst1 = norm * (fabs(w) + fabs(q) + fabs(x)
428                                                                                + fabs(y) + fabs(z));
429                                                                        vr = tst1;
430                                                                        do {
431                                                                                vr = 0.01 * vr;
432                                                                                tst2 = tst1 + vr;
433                                                                        } while (tst2 > tst1);
434                                                                }
435                                                                cdiv(x * r - z * ra + q * sa,
436                                                                        x * s - z * sa - q * ra, vr, vi,
437                                                                     &h[i][na - 1], &h[i][en - 1]);
438                                                                if (fabs(x) > fabs(z) + fabs(q)) {
439                                                                        h[i + 1][na - 1] = (q * h[i][en - 1]
440                                                                                - w * h[i][na - 1] - ra) / x;
441                                                                        h[i + 1][en - 1] = (-sa - w * h[i][en - 1]
442                                                                                - q * h[i][na - 1]) / x;
443                                                                } else
444                                                                        cdiv(-r - y * h[i][na - 1],
445                                                                                -s - y * h[i][en - 1], z, q,
446                                                                             &h[i + 1][na - 1], &h[i + 1][en - 1]);
447                                                        }
448                                                        /* overflow control */
449                                                        t = (fabs(h[i][na-1]) > fabs(h[i][en-1])) ?
450                                                                 fabs(h[i][na-1]) : fabs(h[i][en-1]); /* max */
451                                                        if (t != 0.0) {
452                                                                tst1 = t;
453                                                                tst2 = tst1 + 1.0 / tst1;
454                                                                if (tst2 <= tst1) {
455                                                                        for (j = i; j < en; j++) {
456                                                                                h[j][na - 1] /= t;
457                                                                                h[j][en - 1] /= t;
458                                                                        }
459                                                                }
460                                                        }
461                                                }
462                                        }
463                                }
464                        }
465                }                      /* for nn */
466
467                /* end back substitution. vectors of isolated roots */
468                for (i = 0; i < n; i++) {
469                        if (i + 1 < low || i + 1 > hgh) {
470                                for (j = i; j < n; j++)
471                                        zz[i][j] = h[i][j];
472                        }
473                }
474                /* multiply by transformation matrix to give vectors of
475                 * original full matrix. */
476                for (j = n - 1; j >= low - 1; j--) {
477                        m = ((j + 1) < hgh) ? (j + 1) : hgh; /* min */
478                        for (i = low - 1; i < hgh; i++) {
479                                z = 0.0;
480                                for (k = low - 1; k < m; k++)
481                                        z += zz[i][k] * h[k][j];
482                                zz[i][j] = z;
483                        }
484                }
485        }
486        return;
487
488Lerror: /* two roots found and complex vector */
489        fprintf(stderr, "ERROR in hqr2. two roots found or complex vector");
490        exit(1);
491
492} /*_ hqr2 */
493
494
495void
496readrsrf(r, f, n)
497dmattpmty r;
498dvectpmty f;
499int n;
500{
501        int i, j;
502        double x;
503
504        /* puts("READ 1 Relative Transition Frequencies"); */
505        for (i = 0; i < n; i++) {
506                for (j = 0; j < n; j++) {
507                        fscanf(Rtfifp, "%lf", &x);
508                        if (Debug_optn) printf("%8.1f", x);
509                        r[i][j] = x;
510                }
511                if (Debug_optn) putchar('\n');
512        }
513        if (Debug_optn) putchar('\n');
514        /* puts("READ 2 Relative Transition Frequencies"); */
515        for (i = 0; i < n; i++) {
516                fscanf(Rtfifp, "%lf", &x);
517                if (Debug_optn) printf("%8.1f", x);
518                f[i] = x;
519        }
520        if (Debug_optn) putchar('\n');
521        /* puts("READ # Relative Transition Frequencies"); */
522
523#if 0
524        for (i = 1, x = 0.0; i < Tpmradix; i++) {
525                for (j = 0; j < i; j++) {
526                        printf(" r[%2d][%2d]=%.13e;", i, j, r[i][j]);
527                        if ((j % 2 == 1) || j == (i-1)) putchar('\n');
528                        x += r[i][j];
529                }
530        }
531        for (i = 0; i < Tpmradix; i++) {
532                printf(" f[%2d]=%5.3f;", i, f[i]);
533                if (i % 5 == 4) putchar('\n');
534        }
535        printf("%f\n", x);
536#endif
537
538} /* readrsrf */
539
540
541void
542tpmonepam(a, f)
543dmattpmty a;
544double *f;
545{
546        int i, j, k, tpmhalf;
547        double delta, temp, sum;
548        dvectpmty m;
549
550        for (i = 0, sum = 0.0; i < Tpmradix; i++) {
551                for (j = 0, temp = 0.0; j < Tpmradix; j++)
552                        temp += a[i][j] * f[j];
553                m[i] = temp;
554                sum += temp * f[i];
555        }
556        delta = 0.01 / sum; /* 1 PAM */
557        for (i = 0; i < Tpmradix; i++) {
558                for (j = 0; j < Tpmradix; j++) {
559                        if (i != j)
560                                a[i][j] = delta * a[i][j] * f[j];
561                        else
562#if SH
563                                a[i][j] = - delta * m[i];
564#else
565                                a[i][j] = 1.0 - delta * m[i];
566#endif
567                }
568        }
569
570#ifndef TPM
571        if (Write_optn && !Topting) { /* Write_optn */
572                printf("\nTransition Probability Matrix (x1.0e7)  1PAM\n");
573                if (Tpmradix == NUMAMI) {
574                        tpmhalf = Tpmradix / 2;
575                        for (i = 0; i < Tpmradix; i++) {
576                                printf("%7.0f", a[i][0] * 1.0e7);
577                                for (j = 1; j < tpmhalf; j++)
578                                        printf("%8.0f", a[i][j] * 1.0e7);
579                                putchar('\n');
580                        }
581                        putchar('\n');
582                        for (i = 0; i < Tpmradix; i++) {
583                                printf("%7.0f", a[i][10] * 1.0e7);
584                                for (j = tpmhalf + 1; j < Tpmradix; j++)
585                                        printf("%8.0f", a[i][j] * 1.0e7);
586                                putchar('\n');
587                        }
588                } else {
589                        for (i = 0; i < Tpmradix; i++) {
590                                for (j = 0; j < Tpmradix; j++) {
591                                        if (i == j) {
592                                                for (k = 0, temp = 0; k < Tpmradix; k++) {
593                                                        if (i != k) temp += a[i][k];
594                                                }
595                                                printf("%8.0f", (1.0 - temp) * 1.0e7);
596                                        } else {
597                                                printf("%8.0f", a[i][j]* 1.0e7);
598                                        }
599                                }
600                                putchar('\n');
601                        }
602                }
603        }
604
605        if (Write_optn && Info_optn && !Topting) { /* Write_optn */
606                printf("\nTransition Probability Matrix (x1.0e5)  1PAM\n");
607                printf("%3s", "");
608                for (j = 0; j < Tpmradix; j++) printf("%6s", Cacid3[j]); putchar('\n');
609                for (i = 0; i < Tpmradix; i++) {
610                        printf("%3s", Cacid3[i]);
611                        for (j = 0; j < Tpmradix; j++)
612                                if (i == j) {
613                                        for (k = 0, temp = 0; k < Tpmradix; k++) {
614                                                if (i != k) temp += a[i][k];
615                                        }
616                                        printf("%6.0f", (1.0 - temp) * 1.0e5);
617                                } else {
618                                        printf("%6.0f", a[i][j]* 1.0e5);
619                                }
620                        putchar('\n');
621                }
622                printf("%3s", "Pai");
623                for (j = 0; j < Tpmradix; j++) printf("%6.3f", f[j]); putchar('\n');
624        }
625
626        if (Write_optn && !Topting) { /* Write_optn */
627                if (Tpmradix != NUMAMI) {
628                        /*
629                        printf("\nSubstitution Matrix (x1.0e7)  1PAM\n");
630                        for (i = 0; i < Tpmradix; i++) {
631                                for (j = 0; j < Tpmradix; j++)
632                                        printf("%8.0f", a[i][j] * f[i] * 1.0e7);
633                                putchar('\n');
634                        }
635                        */
636                        printf("\nTS/TV: %.3f\n",
637                                ( a[0][1]*f[0] + a[2][3]*f[2] )
638                                / ( a[0][2]*f[0] + a[0][3]*f[0]
639                                  + a[1][2]*f[1] + a[1][3]*f[1] ) );
640                }
641        }
642#endif /* TPM */
643} /* tpmonepam */
644
645
646void 
647luinverse(omtrx, imtrx, size)
648        dmattpmty omtrx, imtrx;
649        int size;
650{
651        /* INVERSION OF MATRIX ON LU DECOMPOSITION */
652    double eps = 1.0e-20; /* ! */
653        int i, j, k, l, maxi, idx, ix, jx;
654        double sum, tmp, maxb, aw;
655        ivectpmty index;
656        double *wk;
657
658        wk = (double *) malloc((unsigned)size * sizeof(double));
659        aw = 1.0;
660        for (i = 0; i < size; i++) {
661                maxb = 0.0;
662                for (j = 0; j < size; j++) {
663                        if (fabs(omtrx[i][j]) > maxb)
664                                maxb = fabs(omtrx[i][j]);
665                }
666                if (maxb == 0.0) {
667                        fprintf(stderr, "luinverse: singular matrix\n");
668                        exit(1);
669                }
670                wk[i] = 1.0 / maxb;
671        }
672        for (j = 0; j < size; j++) {
673                for (i = 0; i < j; i++) {
674                        sum = omtrx[i][j];
675                        for (k = 0; k < i; k++)
676                                sum -= omtrx[i][k] * omtrx[k][j];
677                        omtrx[i][j] = sum;
678                }
679                maxb = 0.0;
680                for (i = j; i < size; i++) {
681                        sum = omtrx[i][j];
682                        for (k = 0; k < j; k++)
683                                sum -= omtrx[i][k] * omtrx[k][j];
684                        omtrx[i][j] = sum;
685                        tmp = wk[i] * fabs(sum);
686                        if (tmp >= maxb) {
687                                maxb = tmp;
688                                maxi = i;
689                        }
690                }
691                if (j != maxi) {
692                        for (k = 0; k < size; k++) {
693                                tmp = omtrx[maxi][k];
694                                omtrx[maxi][k] = omtrx[j][k];
695                                omtrx[j][k] = tmp;
696                        }
697                        aw = -aw;
698                        wk[maxi] = wk[j];
699                }
700                index[j] = maxi;
701                if (omtrx[j][j] == 0.0)
702                        omtrx[j][j] = eps;
703                if (j != size - 1) {
704                        tmp = 1.0 / omtrx[j][j];
705                        for (i = j + 1; i < size; i++)
706                                omtrx[i][j] *= tmp;
707                }
708        }
709        for (jx = 0; jx < size; jx++) {
710                for (ix = 0; ix < size; ix++)
711                        wk[ix] = 0.0;
712                wk[jx] = 1.0;
713                l = -1;
714                for (i = 0; i < size; i++) {
715                        idx = index[i];
716                        sum = wk[idx];
717                        wk[idx] = wk[i];
718                        if (l != -1) {
719                                for (j = l; j < i; j++)
720                                        sum -= omtrx[i][j] * wk[j];
721                        } else if (sum != 0.0)
722                                l = i;
723                        wk[i] = sum;
724                }
725                for (i = size - 1; i >= 0; i--) {
726                        sum = wk[i];
727                        for (j = i + 1; j < size; j++)
728                                sum -= omtrx[i][j] * wk[j];
729                        wk[i] = sum / omtrx[i][i];
730                }
731                for (ix = 0; ix < size; ix++)
732                        imtrx[ix][jx] = wk[ix];
733        }
734        free((char *)wk);
735        wk = NULL;
736} /*_ luinverse */
737
738
739#ifdef DEBUG
740void 
741mproduct(am, bm, cm, na, nb, nc)
742        dmattpmty am, bm, cm;
743        int na, nb, nc;
744{
745        int ia, ib, ic;
746        double sum;
747
748        for (ia = 0; ia < na; ia++) {
749                for (ic = 0; ic < nc; ic++) {
750                        sum = 0.0;
751                        for (ib = 0; ib < nb; ib++)
752                                sum += am[ia][ib] * bm[ib][ic];
753                        cm[ia][ic] = sum;
754                }
755        }
756} /*_ mproduct */
757
758
759void 
760preigen()
761{
762        int nam1, nam2;
763
764        printf(" EIGEN VECTOR\n");
765        for (nam1 = 0; nam1 < Tpmradix; nam1++) {
766                for (nam2 = 1; nam2 <= Tpmradix; nam2++) {
767                        printf("%7.3f", Evec[nam1][nam2 - 1]);
768                        if (nam2 == 10 || nam2 == 20)
769                                putchar('\n');
770                }
771                putchar('\n');
772        }
773        printf(" EIGEN VALUE\n");
774        for (nam1 = 1; nam1 <= Tpmradix; nam1++) {
775                printf("%7.3f", Eval[nam1 - 1]);
776                if (nam1 == 10 || nam1 == 20)
777                        putchar('\n');
778        }
779} /*_ preigen */
780
781
782void 
783checkevector(imtrx, nn)
784        dmattpmty imtrx;
785        int nn;
786{
787        int i, j;
788
789        for (i = 0; i < nn; i++) {
790                for (j = 0; j < nn; j++) {
791                        if (i == j) {
792                                if (fabs(imtrx[i][j] - 1.0) > 1.0e-5) {
793                                        fprintf(stderr, "ERROR: eigen vector in checkevector 1\n");
794                                        exit(1);
795                                }
796                        } else {
797                                if (fabs(imtrx[i][j]) > 1.0e-5) {
798                                        fprintf(stderr, "ERROR: eigen vector in checkevector 2\n");
799                                        exit(1);
800                                }
801                        }
802                }
803        }
804} /*_ checkevector */
805#endif /* DEBUG */
806
807
808void 
809getrsr(a, ftpm)
810dmattpmty a;
811dvectpmty ftpm;
812{
813        int i, j;
814        double api;
815        dvectpmty forg;
816#ifdef NUC
817        double alpha, beta, beta1, beta2, alpy, alpr;
818#endif /* NUC */
819
820#ifndef NUC
821        api = 0.05;
822#else  /* NUC */
823        api = 0.25;
824#endif /* NUC */
825
826        if (Rrsr_optn) {
827                readrsrf(a, forg, Tpmradix);
828        } else if (Poisn_optn) {
829                for (i = 0; i < Tpmradix; i++) {
830                        for (j = 0; j < Tpmradix; j++) a[i][j] = 1.0;
831                        a[i][i] = 0.0;
832                        forg[i] = api;
833                }
834        } else { /* !Poisn_optn */
835#ifndef NUC
836                if (Mtrev_optn) { /* mtREV */
837                        mtrev(a, forg);
838                } else if (Dayhf_optn) { /* Dayhoff */
839                        dyhfjtt(a, forg, TRUE);
840                } else { /* JTT */
841                        dyhfjtt(a, forg, FALSE);
842                }
843#else   /* NUC */
844                beta = 1.0;
845                beta1 = (beta * 2.0) / (Beta12 + 1.0);
846                beta2 = Beta12 * beta1;
847                alpha = AlphaBeta;
848                alpr  = (alpha * 2.0) / (AlphaYR + 1.0);
849                alpy  = AlphaYR * alpr;
850        a[0][0] =  0.0;  a[0][1] = alpy;  a[0][2] = beta1; a[0][3] = beta1;
851                a[1][0] = alpy;  a[1][1] =  0.0;  a[1][2] = beta2; a[1][3] = beta2;
852                a[2][0] = beta1; a[2][1] = beta2; a[2][2] =  0.0;  a[2][3] = alpr;
853                a[3][0] = beta1; a[3][1] = beta2; a[3][2] = alpr;  a[3][3] =  0.0;
854                forg[0] = 0.25;  forg[1] = 0.25;  forg[2] = 0.25;  forg[3] = 0.25;
855#endif  /* NUC */
856        } /* Poisn_optn */
857
858#ifdef DEBUG
859        for (i = 0; i < Tpmradix; i++) {
860                for (j = 0; j < Tpmradix; j++) {
861                        if (a[i][j] != a[j][i]) {
862                                fputs("ERROR: Relative Substitution Rate Matrix: ", stderr);
863                                fprintf(stderr, "a[%d][%d] = %.3f != a[%d][%d]\n",
864                                        i, j, a[i][j], j, i);
865                                exit(1);
866                        }
867                }
868        }
869#endif /* DEBUG */
870
871        if (Write_optn && !Topting) { /* Debug_optn */
872                puts("\nRelative Substitution Rate Matrix");
873                for (i = 0; i < Tpmradix; i++) {
874                        for (j = 0; j < Tpmradix; j++)
875#ifndef NUC
876                                if (i != j)
877                                        printf("%5.0f", a[i][j]);
878                                else
879                                        printf("%5s", Cacid3[i]);
880#else                   /* NUC */
881                                if (i != j)
882                                        printf("%8.3f", a[i][j]);
883                                else
884                                        printf("%6s  ", Cacid1[i]);
885#endif                  /* NUC */
886                        putchar('\n');
887                }
888        }
889
890#ifdef TPM
891        for (i = 0; i < Tpmradix; i++) {
892                for (j = 0; j < Tpmradix; j++)
893                        Rtf[i][j] = a[i][j];
894        }
895#endif /* TPM */
896
897        if (Frequ_optn) {
898                for (i = 0; i < Tpmradix - 1; i++) {
899                        for (j = i + 1; j < Tpmradix; j++) {
900                                if (Freqemp[i] == Freqemp[j]) {
901                                        Freqemp[i] += 0.00001;
902                                        Freqemp[j] -= 0.00001;
903                                }
904                        }
905                }
906                for (i = 0; i < Tpmradix; i++) ftpm[i] = Freqemp[i];
907        } else {
908                for (i = 0; i < Tpmradix; i++) ftpm[i] = forg[i];
909        }
910} /* getrsr */
911
912
913void 
914tranprobmat()
915{
916        /* make transition probability matrix */
917        dmattpmty a, b;
918        dvectpmty evali;
919        ivectpmty ordr;
920        int i, j, k, err;
921        double zero, temp;
922        boolean error;
923
924        getrsr(a, Freqtpm);
925
926        tpmonepam(a, Freqtpm);
927
928        for (i = 0; i < Tpmradix; i++) {
929                for (j = 0; j < Tpmradix; j++)
930                        b[i][j] = a[i][j];
931        }
932        error = FALSE;
933        elmhes(a, ordr, Tpmradix);
934        eltran(a, Evec, ordr, Tpmradix);
935        hqr2(Tpmradix, 1, Tpmradix, &err, a, Evec, Eval, evali);
936
937#ifdef DEBUG
938        for (j = 0; j < Tpmradix; j++) {
939                for (i = 0, zero = 0.0; i < Tpmradix; i++) {
940                        for (k = 0; k < Tpmradix; k++)
941                                zero += b[i][k] * Evec[k][j];
942                        zero -= Eval[j] * Evec[i][j];
943                        if (fabs(zero) > 1.0e-5) {
944                                error = TRUE;
945                                printf(" ERROR: Eigen System%  .3E", zero);
946                                if (i % 5 == 0) putchar('\n');
947                        }
948                }
949        }
950        if (error) {
951                fflush(stdout);
952                fputs("\nERROR: Eigen Value of Transition Probability Matrix\n",stderr);
953                exit(1);
954        }
955#endif /* DEBUG */
956
957        for (i = 0; i < Tpmradix; i++) {
958                for (j = 0; j < Tpmradix; j++)
959                        b[i][j] = Evec[i][j];
960        }
961        luinverse(b, Ievc, Tpmradix);
962
963#ifdef DEBUG
964        mproduct(Evec, Ievc, b, Tpmradix, Tpmradix, Tpmradix);
965        checkevector(b, Tpmradix);
966        if (Debug_optn) preigen();
967#endif /* DEBUG */
968
969        for (i = 0; i < Tpmradix - 1; i++) { /* !? */
970                for (j = i + 1; j < Tpmradix; j++) {
971                        temp       = Ievc[i][j];
972                        Ievc[i][j] = Ievc[j][i];
973                        Ievc[j][i] = temp;
974                }
975        }
976#if SH
977        for (i = 0; i < Tpmradix; i++) Evl2[i] = Eval[i] * Eval[i];
978#else
979        for (i = 0; i < Tpmradix; i++) {
980                temp = log(Eval[i]);
981                Eval[i] = temp;
982                Evl2[i] = temp * temp;
983        }
984#endif
985} /*_ tranprobmat */
986
987
988#if VARITPM
989void 
990varitpm()
991{
992        /* make transition probability matrix */
993        dmattpmty a, b, r, o, rr;
994        dvectpmty evali;
995        ivectpmty ordr;
996        int i, j, k, err, ii, jj, kk;
997        double zero, temp;
998        double aij, lklorg, lklnew, d1, d2, dd, vari, notch;
999        double tlo, tro, tln, trn, tl1, tl2, tr1, tr2, tld, trd, ttl, ttr, varil, varir;
1000
1001        notch = Percent;
1002        initpartlkl(Ctree);
1003        aproxlkl(Ctree);
1004        lklorg = Ctree->aproxl;
1005
1006        getrsr(r, Freqtpm);
1007        for (i = 0; i < Tpmradix; i++) {
1008                for (j = 0; j < Tpmradix; j++) o[i][j] = r[i][j];
1009        }
1010        tpmonepam(o, Freqtpm);
1011        for (j = 0; j < Tpmradix; j++) rr[j][j] = 0.0;
1012
1013        for (ii = 1; ii < Tpmradix; ii++) {
1014        for (jj = 0; jj < ii; jj++) {
1015
1016        aij = r[ii][jj];
1017        tlo = o[ii][jj];
1018        tro = o[jj][ii];
1019
1020        for (kk = 0; kk < 2; kk++) {
1021
1022                for (i = 0; i < Tpmradix; i++) {
1023                        for (j = 0; j < Tpmradix; j++) a[i][j] = r[i][j];
1024                }
1025       
1026                if (kk == 0) {
1027                        a[ii][jj] -= notch;
1028                        a[jj][ii] -= notch;
1029                } else {
1030                        a[ii][jj] += notch;
1031                        a[jj][ii] += notch;
1032                }
1033       
1034                tpmonepam(a, Freqtpm);
1035       
1036                tln = a[ii][jj];
1037                trn = a[jj][ii];
1038       
1039                for (i = 0; i < Tpmradix; i++) {
1040                        for (j = 0; j < Tpmradix; j++)
1041                                b[i][j] = a[i][j];
1042                }
1043                elmhes(a, ordr, Tpmradix);
1044                eltran(a, Evec, ordr, Tpmradix);
1045                hqr2(Tpmradix, 1, Tpmradix, &err, a, Evec, Eval, evali);
1046       
1047                for (i = 0; i < Tpmradix; i++) {
1048                        for (j = 0; j < Tpmradix; j++)
1049                                b[i][j] = Evec[i][j];
1050                }
1051                luinverse(b, Ievc, Tpmradix);
1052       
1053                for (i = 0; i < Tpmradix - 1; i++) { /* !? */
1054                        for (j = i + 1; j < Tpmradix; j++) {
1055                                temp       = Ievc[i][j];
1056                                Ievc[i][j] = Ievc[j][i];
1057                                Ievc[j][i] = temp;
1058                        }
1059                }
1060#if SH
1061                for (i = 0; i < Tpmradix; i++) Evl2[i] = Eval[i] * Eval[i];
1062#else
1063                for (i = 0; i < Tpmradix; i++) {
1064                        temp = log(Eval[i]);
1065                        Eval[i] = temp;
1066                        Evl2[i] = temp * temp;
1067                }
1068#endif
1069       
1070                initpartlkl(Ctree);
1071                aproxlkl(Ctree);
1072                lklnew = Ctree->aproxl;
1073                if (kk == 0) {
1074                        d1  = (lklorg - lklnew) / notch;
1075                        tl1 = (lklorg - lklnew) / (tlo - tln); 
1076                        tr1 = (lklorg - lklnew) / (tro - trn);
1077                        tld = (tlo - tln) / 2.0;
1078                        trd = (tro - trn) / 2.0;
1079                } else {
1080                        d2  = (lklnew - lklorg) / notch;
1081                        tl2 = (lklorg - lklnew) / (tln - tlo); 
1082                        tr2 = (lklorg - lklnew) / (trn - tro);
1083                        tld += (tln - tlo) / 2.0;
1084                        trd += (trn - tro) / 2.0;
1085                }
1086                /*
1087                printf("%9.3f%9.3f%9.3f%9.3f%9.6f%9.6f\n",
1088                        tlo*1e5, tln*1e5, tro*1e5, trn*1e5, tld*1e5, trd*1e5);
1089                */
1090
1091        }
1092        dd  = (d2 - d1) / notch;
1093        ttl = (tl2 - tl1) / tld;
1094        ttr = (tr2 - tr1) / trd;
1095        vari  =  sqrt(1.0 / fabs(dd));
1096        varil =  sqrt(1.0 / fabs(ttl));
1097        varir =  sqrt(1.0 / fabs(ttr));
1098        printf("%3d%3d%2s%2s%9.5f%9.5f%3s%6.0f%6.0f %8.3f%8.3f%8.3f%8.3f\n",
1099                ii+1, jj+1, Cacid1[ii], Cacid1[jj], d1, d2, (dd < 0.0 ? "mi" : "PL"),
1100                aij, vari, tlo*100000, varil*100000, tro*100000, varir*100000);
1101        rr[ii][jj] = aij;
1102        if (dd < 0.0) {
1103                rr[jj][ii] = vari;
1104        } else {
1105                rr[jj][ii] = 0.0;
1106        }
1107        }
1108        }
1109        printf("\nnotch:%12.8f\n", notch);
1110        for (i = 0; i < Tpmradix; i++) {
1111                for (j = 0; j < Tpmradix; j++) {
1112                        if (j == i) {
1113                                printf("%5.3s", Cacid3[i]);
1114                        } else {
1115                                if (i > j) {
1116                                        printf("%5.0f", rr[i][j]);
1117                                } else {
1118                                        if (rr[j][i] == 1.0)
1119                                                printf("%5.1s", "-");
1120                                        else if (rr[i][j] == 0.0)
1121                                                printf("%5.1s", "*");
1122                                        else
1123                                                printf("%5.0f", rr[i][j]);
1124                                }
1125                        }
1126                } putchar('\n');
1127        } 
1128} /*_ varitpm */
1129#endif /* VARITPM */
1130
1131
1132void
1133prfreq()
1134{
1135        int i;
1136
1137        printf("\n%s Acid Frequencies\n", Infomol);
1138        printf("%3s%4s%8s%7s%7s\n",
1139                "", "", "Model ", "Data", "     ");
1140        for (i = 0; i < Tpmradix; i++) {
1141                printf("%3d%3s%8.3f%8.3f\n",
1142                        i + 1, Cacid1[i], Freqtpm[i], Freqemp[i]);
1143        }
1144} /*_ refreq */
1145
1146
1147void 
1148tprobmtrx(arc, tpr)
1149        double arc;
1150        dmattpmty tpr;
1151{
1152        /* TRANSITION PROBABILITY MATRIX */
1153        register int i, j, k;
1154        register double temp;
1155        dvectpmty vexp;
1156        dmattpmty iexp;
1157
1158        for (k = 0; k < Tpmradix; k++)
1159                vexp[k] = exp(arc * Eval[k]);
1160        for (j = 0; j < Tpmradix; j++) {
1161                for (k = 0; k < Tpmradix; k++)
1162                        iexp[j][k] = Ievc[j][k] * vexp[k];
1163        }
1164
1165        for (i = 0; i < Tpmradix; i++) {
1166                for (j = 0; j < Tpmradix; j++) {
1167                        temp = 0.0;
1168                        for (k = 0; k < Tpmradix; k++)
1169                                temp += Evec[i][k] * iexp[j][k];
1170                        tpr[i][j] = temp;
1171                }
1172        }
1173} /*_ tprobmtrx */
1174
1175void 
1176tprobmtrxt(arc, tpr)
1177        double arc;
1178        dmattpmty tpr;
1179{
1180        /* TRANSITION PROBABILITY MATRIX */
1181        register int i, j, k;
1182        register double temp;
1183        dvectpmty vexp;
1184        dmattpmty iexp;
1185
1186        for (k = 0; k < Tpmradix; k++)
1187                vexp[k] = exp(arc * Eval[k]);
1188        for (j = 0; j < Tpmradix; j++) {
1189                for (k = 0; k < Tpmradix; k++)
1190                        iexp[j][k] = Ievc[j][k] * vexp[k];
1191        }
1192
1193        for (i = 0; i < Tpmradix; i++) {
1194                for (j = 0; j < Tpmradix; j++) {
1195                        temp = 0.0;
1196                        for (k = 0; k < Tpmradix; k++)
1197                                temp += Evec[i][k] * iexp[j][k];
1198                        tpr[j][i] = temp; /* i <-> j */
1199                }
1200        }
1201} /*_ tprobmtrxt */
1202
1203
1204void 
1205tdiffmtrx(arc, tpr, td1, td2)
1206        double arc;
1207        dmattpmty tpr, td1, td2;
1208{
1209        /* TRANSITION PROBABILITY AND DIFFRENCE MATRIX */
1210        register int i, j, k;
1211        register double x, tp, t1, t2;
1212        dvectpmty vexp;
1213        dmattpmty iexp;
1214        dvector evc;
1215
1216        for (k = 0; k < Tpmradix; k++)
1217                vexp[k] = exp(arc * Eval[k]);
1218        for (j = 0; j < Tpmradix; j++) {
1219                for (k = 0; k < Tpmradix; k++)
1220                        iexp[j][k] = Ievc[j][k] * vexp[k];
1221        }
1222
1223        for (i = 0; i < Tpmradix; i++) {
1224                evc = Evec[i];
1225                for (j = 0; j < Tpmradix; j++) {
1226                        tp = t1 = t2 = 0.0;
1227                        for (k = 0; k < Tpmradix; k++) {
1228                                x = evc[k] * iexp[j][k];
1229                                tp += x;
1230                                t1 += x * Eval[k];
1231                                t2 += x * Evl2[k];
1232                        }
1233                        tpr[i][j] = tp;
1234                        td1[i][j] = t1;
1235                        td2[i][j] = t2;
1236                }
1237        }
1238} /*_ tdiffmtrx */
1239
1240
1241#if 0
1242void
1243tprobmtrx2(arc, tpr) /* Ievc is not Transposition matrix */
1244        double arc;
1245        dmattpmty tpr;
1246{
1247        /* TRANSITION PROBABILITY MATRIX */
1248        register int i, j, k;
1249        register double s0, s1, s2, s3, s4, s5, s6, s7;
1250        register double s8, s9, s10, s11, s12, s13, s14, s15;
1251        register double b0, b1, b2, b3, c0, c1, c2, c3, v0, v1;
1252        dmattpmty iexp;
1253
1254        for (k = 0; k < Tpmradix; k++) {
1255                s0 = exp(arc * Eval[k]);
1256                for (j = 0; j < Tpmradix; j++)
1257                        iexp[k][j] = Ievc[k][j] * s0;
1258        }
1259        for (i = 0; i < Tpmradix; i += 4) {
1260                for (j = 0; j < Tpmradix; j += 4) {
1261                        s0  = 0.0; s1  = 0.0; s2  = 0.0; s3  = 0.0;
1262                        s4  = 0.0; s5  = 0.0; s6  = 0.0; s7  = 0.0;
1263                        s8  = 0.0; s9  = 0.0; s10 = 0.0; s11 = 0.0;
1264                        s12 = 0.0; s13 = 0.0; s14 = 0.0; s15 = 0.0;
1265                        b0 = Evec[i  ][0];
1266                        b1 = Evec[i+1][0];
1267                        c0 = iexp[0][j];
1268                        c1 = iexp[0][j+1];
1269                        c2 = iexp[0][j+2];
1270                        c3 = iexp[0][j+3];
1271                        v0 = b0 * c0;
1272                        v1 = b0 * c1;
1273                        for (k = 0; k < Tpmradix-1; k++) {
1274                                s0  += v0;
1275                                s4  += v1;
1276                                s8  += b0 * c2;
1277                                s12 += b0 * c3;
1278                                s1  += b1 * c0;
1279                                s5  += b1 * c1;
1280                                s9  += b1 * c2;
1281                                s13 += b1 * c3;
1282                                b0 = Evec[i  ][k+1];
1283                                b1 = Evec[i+1][k+1];
1284                                b2 = Evec[i+2][k];
1285                                b3 = Evec[i+3][k];
1286                                s2  += b2 * c0;
1287                                s6  += b2 * c1;
1288                                s10 += b2 * c2;
1289                                s14 += b2 * c3;
1290                                s3  += b3 * c0;
1291                                s7  += b3 * c1;
1292                                s11 += b3 * c2;
1293                                s15 += b3 * c3;
1294                                c0 = iexp[k+1][j];
1295                                c1 = iexp[k+1][j+1];
1296                                c2 = iexp[k+1][j+2];
1297                                c3 = iexp[k+1][j+3];
1298                                v0 = b0 * c0;
1299                                v1 = b0 * c1;
1300                        }
1301                        s0  += v0;
1302                        s4  += v1;
1303                        s8  += b0 * c2;
1304                        s12 += b0 * c3;
1305                        s1  += b1 * c0;
1306                        s5  += b1 * c1;
1307                        s9  += b1 * c2;
1308                        s13 += b1 * c3;
1309                        b2 = Evec[i+2][k];
1310                        b3 = Evec[i+3][k];
1311                        s2  += b2 * c0;
1312                        s6  += b2 * c1;
1313                        s10 += b2 * c2;
1314                        s14 += b2 * c3;
1315                        s3  += b3 * c0;
1316                        s7  += b3 * c1;
1317                        s11 += b3 * c2;
1318                        s15 += b3 * c3;
1319                        tpr[i  ][j]   = s0;
1320                        tpr[i  ][j+1] = s4;
1321                        tpr[i  ][j+2] = s8;
1322                        tpr[i  ][j+3] = s12;
1323                        tpr[i+1][j]   = s1;
1324                        tpr[i+1][j+1] = s5;
1325                        tpr[i+1][j+2] = s9;
1326                        tpr[i+1][j+3] = s13;
1327                        tpr[i+2][j]   = s2;
1328                        tpr[i+2][j+1] = s6;
1329                        tpr[i+2][j+2] = s10;
1330                        tpr[i+2][j+3] = s14;
1331                        tpr[i+3][j]   = s3;
1332                        tpr[i+3][j+1] = s7;
1333                        tpr[i+3][j+2] = s11;
1334                        tpr[i+3][j+3] = s15;
1335                }
1336        }
1337} /*_ tprobmtrx */
1338#endif
1339
1340
1341#if 0
1342void
1343tdiffmtrx2(arc, tpr, td1, td2) /* Ievc is not Transposition matrix */
1344        double arc;
1345        dmattpmty tpr, td1, td2;
1346{
1347        /* TRANSITION PROBABILITY MATRIX */
1348        register int i, j, k;
1349        register double s00, s01, s10, s11;
1350        register double t00, t01, t10, t11;
1351        register double u00, u01, u10, u11;
1352        register double v00, v01, v10, v11;
1353        register double b0, b1, c0, c1;
1354        register double x, y, z;
1355        dvectpmty ex0, ex1, ex2;
1356
1357        for (k = 0; k < Tpmradix; k++) {
1358                x = exp(arc * Eval[k]);
1359                ex0[k] = x;
1360                ex1[k] = Eval[k] * x;
1361                ex2[k] = Evl2[k] * x;
1362        }
1363
1364        for (i = 0; i < Tpmradix; i += 2) {
1365                for (j = 0; j < Tpmradix; j += 2) {
1366                        s00 = 0.0; s01 = 0.0; s10 = 0.0; s11 = 0.0;
1367                        t00 = 0.0; t01 = 0.0; t10 = 0.0; t11 = 0.0;
1368                        u00 = 0.0; u01 = 0.0; u10 = 0.0; u11 = 0.0;
1369
1370                        b0 = Evec[i][0];
1371                        c0 = Ievc[0][j];
1372                        c1 = Ievc[0][j+1];
1373                        for (k = 0; k < Tpmradix-1; k++) {
1374                                x = ex0[k];
1375                                y = ex1[k];
1376                                z = ex2[k];
1377                                v00 = b0 * c0;
1378                                v01 = b0 * c1;
1379                                s00 += v00 * x;
1380                                s01 += v01 * x;
1381                                t00 += v00 * y;
1382                                t01 += v01 * y;
1383                                u00 += v00 * z;
1384                                u01 += v01 * z;
1385                                b0 = Evec[i][k+1];
1386                                b1 = Evec[i+1][k];
1387                                v10 = b1 * c0;
1388                                v11 = b1 * c1;
1389                                s10 += v10 * x;
1390                                s11 += v11 * x;
1391                                t10 += v10 * y;
1392                                t11 += v11 * y;
1393                                u10 += v10 * z;
1394                                u11 += v11 * z;
1395                                c0 = Ievc[k+1][j];
1396                                c1 = Ievc[k+1][j+1];
1397                        }
1398                        x = ex0[Tpmradix-1];
1399                        y = ex1[Tpmradix-1];
1400                        z = ex2[Tpmradix-1];
1401                                v00 = b0 * c0;
1402                                v01 = b0 * c1;
1403                                s00 += v00 * x;
1404                                s01 += v01 * x;
1405                                t00 += v00 * y;
1406                                t01 += v01 * y;
1407                                u00 += v00 * z;
1408                                u01 += v01 * z;
1409                                b1 = Evec[i+1][k];
1410                                v10 = b1 * c0;
1411                                v11 = b1 * c1;
1412                                s10 += v10 * x;
1413                                s11 += v11 * x;
1414                                t10 += v10 * y;
1415                                t11 += v11 * y;
1416                                u10 += v10 * z;
1417                                u11 += v11 * z;
1418                        tpr[i  ][j  ] = s00;
1419                        tpr[i  ][j+1] = s01;
1420                        tpr[i+1][j  ] = s10;
1421                        tpr[i+1][j+1] = s11;
1422                        td1[i  ][j  ] = t00;
1423                        td1[i  ][j+1] = t01;
1424                        td1[i+1][j  ] = t10;
1425                        td1[i+1][j+1] = t11;
1426                        td2[i  ][j  ] = u00;
1427                        td2[i  ][j+1] = u01;
1428                        td2[i+1][j  ] = u10;
1429                        td2[i+1][j+1] = u11;
1430                }
1431        }
1432} /*_ tdiffmtrx */
1433#endif
Note: See TracBrowser for help on using the repository browser.