dopla.c

完全使用C++写的高效线性代数运算库！还提供了矩阵类。

/*        Note:  The operation count is overestimated because */
/*        it is assumed that the factor U fills in to the maximum */
/*        extent, i.e., that its bandwidth goes from KU to KL + KU. */

    } else if (lsamen_(&c__2, c2, "GB", 2L, 2L)) {

/*        xGBTRF:  M, N, KL, KU  =>  M, N, KL, KU */

    if (lsamen_(&c__3, c3, "TRF", 3L, 3L)) {
        for (i = min(*m,*n); i >= 1; --i) {
/* Computing MAX */
/* Computing MIN */
        i__3 = *kl, i__4 = *m - i;
        i__1 = 0, i__2 = min(i__3,i__4);
        wl = (doublereal) max(i__1,i__2);
/* Computing MAX */
/* Computing MIN */
        i__3 = *kl + *ku, i__4 = *n - i;
        i__1 = 0, i__2 = min(i__3,i__4);
        wu = (doublereal) max(i__1,i__2);
        mults += wl * (wu + 1.);
        adds += wl * wu;
/* L10: */
        }

/*        xGBTRS:  N, NRHS, KL, KU  =>  M, N, KL, KU */

    } else if (lsamen_(&c__3, c3, "TRS", 3L, 3L)) {
/* Computing MAX */
/* Computing MIN */
        i__3 = *kl, i__4 = *m - 1;
        i__1 = 0, i__2 = min(i__3,i__4);
        wl = (doublereal) max(i__1,i__2);
/* Computing MAX */
/* Computing MIN */
        i__3 = *kl + *ku, i__4 = *m - 1;
        i__1 = 0, i__2 = min(i__3,i__4);
        wu = (doublereal) max(i__1,i__2);
        mults = en * (em * (wl + 1. + wu) - (wl * (wl + 1.) + wu * (wu + 
            1.)) * .5);
        adds = en * (em * (wl + wu) - (wl * (wl + 1.) + wu * (wu + 1.)) * 
            .5);

    }

/*     -------------------------------------- */
/*     PO:  POsitive definite matrices */
/*     PP:  Positive definite Packed matrices */
/*     -------------------------------------- */

    } else if (lsamen_(&c__2, c2, "PO", 2L, 2L) || lsamen_(&c__2, c2, "PP", 
        2L, 2L)) {

/*        xPOTRF:  N  =>  M */

    if (lsamen_(&c__3, c3, "TRF", 3L, 3L)) {
        mults = em * (em * (em * .16666666666666666 + .5) + 
            .33333333333333331);
        adds = em * .16666666666666666 * (em * em - 1.);

/*        xPOTRS:  N, NRHS  =>  M, N */

    } else if (lsamen_(&c__3, c3, "TRS", 3L, 3L)) {
        mults = en * (em * (em + 1.));
        adds = en * (em * (em - 1.));

/*        xPOTRI:  N  =>  M */

    } else if (lsamen_(&c__3, c3, "TRI", 3L, 3L)) {
        mults = em * (em * (em * .33333333333333331 + 1.) + 
            .66666666666666663);
        adds = em * (em * (em * .33333333333333331 - .5) + 
            .16666666666666666);

    }

/*     ------------------------------------ */
/*     PB:  Positive definite Band matrices */
/*     ------------------------------------ */

    } else if (lsamen_(&c__2, c2, "PB", 2L, 2L)) {

/*        xPBTRF:  N, K  =>  M, KL */

    if (lsamen_(&c__3, c3, "TRF", 3L, 3L)) {
        mults = ek * (ek * (ek * -.33333333333333331 - 1.) - 
            .66666666666666663) + em * (ek * (ek * .5 + 1.5) + 1.);
        adds = ek * (ek * (ek * -.33333333333333331 - .5) - 
            .16666666666666666) + em * (ek / 2. * (ek + 1.));

/*        xPBTRS:  N, NRHS, K  =>  M, N, KL */

    } else if (lsamen_(&c__3, c3, "TRS", 3L, 3L)) {
        mults = en * ((em * 2 - ek) * (ek + 1.));
        adds = en * (ek * (em * 2 - (ek + 1.)));

    }

/*     -------------------------------------------------------- */
/*     SY:  SYmmetric indefinite matrices */
/*     SP:  Symmetric indefinite Packed matrices */
/*     HE:  HErmitian indefinite matrices (complex only) */
/*     HP:  Hermitian indefinite Packed matrices (complex only) */
/*     -------------------------------------------------------- */

    } else if (lsamen_(&c__2, c2, "SY", 2L, 2L) || lsamen_(&c__2, c2, "SP", 
        2L, 2L) || lsamen_(&c__3, subnam, "ZHE", 6L, 3L) || lsamen_(&c__3,
         subnam, "ZHE", 6L, 3L) || lsamen_(&c__3, subnam, "ZHP", 6L, 3L) 
        || lsamen_(&c__3, subnam, "ZHP", 6L, 3L)) {

/*        xSYTRF:  N  =>  M */

    if (lsamen_(&c__3, c3, "TRF", 3L, 3L)) {
        mults = em * (em * (em * .16666666666666666 + .5) + 
            3.3333333333333335);
        adds = em / 6. * (em * em - 1.);

/*        xSYTRS:  N, NRHS  =>  M, N */

    } else if (lsamen_(&c__3, c3, "TRS", 3L, 3L)) {
        mults = en * em * em;
        adds = en * (em * (em - 1.));

/*        xSYTRI:  N  =>  M */

    } else if (lsamen_(&c__3, c3, "TRI", 3L, 3L)) {
        mults = em * (em * em * .33333333333333331 + .66666666666666663);
        adds = em * (em * em * .33333333333333331 - .33333333333333331);

/*        xSYTRD, xSYTD2:  N  =>  M */

    } else if (lsamen_(&c__3, c3, "TRD", 3L, 3L) || lsamen_(&c__3, c3, 
        "TD2", 3L, 3L)) {
        if (*m == 1) {
        mults = 0.;
        adds = 0.;
        } else {
        mults = em * (em * (em * .66666666666666663 + 2.5) - 
            .16666666666666666) - 15.;
        adds = em * (em * (em * .66666666666666663 + 1.) - 
            2.6666666666666665) - 4.;
        }
    }

/*     ------------------- */
/*     Triangular matrices */
/*     ------------------- */

    } else if (lsamen_(&c__2, c2, "TR", 2L, 2L) || lsamen_(&c__2, c2, "TP", 
        2L, 2L)) {

/*        xTRTRS:  N, NRHS  =>  M, N */

    if (lsamen_(&c__3, c3, "TRS", 3L, 3L)) {
        mults = en * em * (em + 1.) / 2.;
        adds = en * em * (em - 1.) / 2.;

/*        xTRTRI:  N  =>  M */

    } else if (lsamen_(&c__3, c3, "TRI", 3L, 3L)) {
        mults = em * (em * (em * .16666666666666666 + .5) + 
            .33333333333333331);
        adds = em * (em * (em * .16666666666666666 - .5) + 
            .33333333333333331);

    }

    } else if (lsamen_(&c__2, c2, "TB", 2L, 2L)) {

/*        xTBTRS:  N, NRHS, K  =>  M, N, KL */

    if (lsamen_(&c__3, c3, "TRS", 3L, 3L)) {
        mults = en * (em * (em + 1.) / 2. - (em - ek - 1.) * (em - ek) / 
            2.);
        adds = en * (em * (em - 1.) / 2. - (em - ek - 1.) * (em - ek) / 
            2.);
    }

/*     -------------------- */
/*     Trapezoidal matrices */
/*     -------------------- */

    } else if (lsamen_(&c__2, c2, "TZ", 2L, 2L)) {

/*        xTZRQF:  M, N => M, N */

    if (lsamen_(&c__3, c3, "RQF", 3L, 3L)) {
        emn = (doublereal) min(*m,*n);
        mults = em * 3 * (en - em + 1) + (en * 2 - em * 2 + 3) * (em * em 
            - emn * (emn + 1) / 2);
        adds = (en - em + 1) * (em + em * 2 * em - emn * (emn + 1));
    }

/*     ------------------- */
/*     Orthogonal matrices */
/*     ------------------- */

    } else if (sord && lsamen_(&c__2, c2, "OR", 2L, 2L) || corz && lsamen_(&
        c__2, c2, "UN", 2L, 2L)) {

/*        -MQR, -MLQ, -MQL, or -MRQ:  M, N, K, SIDE  =>  M, N, KL, KU 
*/
/*           where KU<= 0 indicates SIDE = 'L' */
/*           and   KU> 0  indicates SIDE = 'R' */

    if (lsamen_(&c__3, c3, "MQR", 3L, 3L) || lsamen_(&c__3, c3, "MLQ", 3L,
         3L) || lsamen_(&c__3, c3, "MQL", 3L, 3L) || lsamen_(&c__3, 
        c3, "MRQ", 3L, 3L)) {
        if (*ku <= 0) {
        mults = ek * en * (em * 2. + 2. - ek);
        adds = ek * en * (em * 2. + 1. - ek);
        } else {
        mults = ek * (em * (en * 2. - ek) + (em + en + (1. - ek) / 2.)
            );
        adds = ek * em * (en * 2. + 1. - ek);
        }

/*        -GQR or -GQL:  M, N, K  =>  M, N, KL */

    } else if (lsamen_(&c__3, c3, "GQR", 3L, 3L) || lsamen_(&c__3, c3, 
        "GQL", 3L, 3L)) {
        mults = ek * (en * 2. - ek - 1.6666666666666667 + (em * 2. * en + 
            ek * (ek * .66666666666666663 - em - en)));
        adds = ek * (en - em + .33333333333333331 + (em * 2. * en + ek * (
            ek * .66666666666666663 - em - en)));

/*        -GLQ or -GRQ:  M, N, K  =>  M, N, KL */

    } else if (lsamen_(&c__3, c3, "GLQ", 3L, 3L) || lsamen_(&c__3, c3, 
        "GRQ", 3L, 3L)) {
        mults = ek * (em + en - ek - .66666666666666663 + (em * 2. * en + 
            ek * (ek * .66666666666666663 - em - en)));
        adds = ek * (em - en + .33333333333333331 + (em * 2. * en + ek * (
            ek * .66666666666666663 - em - en)));

    }

    }

    ret_val = mulfac * mults + addfac * adds;

    return ret_val;

/*     End of DOPLA */

} /* dopla_ */