dopla.c

提供矩阵类的函数库

#include "blaswrap.h"
/*  -- translated by f2c (version 19990503).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Table of constant values */

static integer c__2 = 2;
static integer c__3 = 3;

doublereal dopla_(char *subnam, integer *m, integer *n, integer *kl, integer *
	ku, integer *nb)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;
    doublereal ret_val;

    /* Builtin functions   
       Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);

    /* Local variables */
    static doublereal adds;
    static logical sord, corz;
    static integer i__;
    extern logical lsame_(char *, char *);
    static char c1[1], c2[2], c3[3];
    static doublereal mults, addfac, ek, em, en, wl, mulfac, wu;
    extern logical lsamen_(integer *, char *, char *);
    static doublereal emn;


/*  -- LAPACK timing routine (version 3.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       June 30, 1999   


    Purpose   
    =======   

    DOPLA computes an approximation of the number of floating point   
    operations used by the subroutine SUBNAM with the given values   
    of the parameters M, N, KL, KU, and NB.   

    This version counts operations for the LAPACK subroutines.   

    Arguments   
    =========   

    SUBNAM  (input) CHARACTER*6   
            The name of the subroutine.   

    M       (input) INTEGER   
            The number of rows of the coefficient matrix.  M >= 0.   

    N       (input) INTEGER   
            The number of columns of the coefficient matrix.   
            For solve routine when the matrix is square,   
            N is the number of right hand sides.  N >= 0.   

    KL      (input) INTEGER   
            The lower band width of the coefficient matrix.   
            If needed, 0 <= KL <= M-1.   
            For xGEQRS, KL is the number of right hand sides.   

    KU      (input) INTEGER   
            The upper band width of the coefficient matrix.   
            If needed, 0 <= KU <= N-1.   

    NB      (input) INTEGER   
            The block size.  If needed, NB >= 1.   

    Notes   
    =====   

    In the comments below, the association is given between arguments   
    in the requested subroutine and local arguments.  For example,   

    xGETRS:  N, NRHS  =>  M, N   

    means that arguments N and NRHS in DGETRS are passed to arguments   
    M and N in this procedure.   

    =====================================================================   


       --------------------------------------------------------   
       Initialize DOPLA to 0 and do a quick return if possible.   
       -------------------------------------------------------- */

    ret_val = 0.;
    mults = 0.;
    adds = 0.;
    *(unsigned char *)c1 = *(unsigned char *)subnam;
    s_copy(c2, subnam + 1, (ftnlen)2, (ftnlen)2);
    s_copy(c3, subnam + 3, (ftnlen)3, (ftnlen)3);
    sord = lsame_(c1, "S") || lsame_(c1, "D");
    corz = lsame_(c1, "C") || lsame_(c1, "Z");
    if (*m <= 0 || ! (sord || corz)) {
	return ret_val;
    }

/*     ---------------------------------------------------------   
       If the coefficient matrix is real, count each add as 1   
       operation and each multiply as 1 operation.   
       If the coefficient matrix is complex, count each add as 2   
       operations and each multiply as 6 operations.   
       --------------------------------------------------------- */

    if (lsame_(c1, "S") || lsame_(c1, "D")) {
	addfac = 1.;
	mulfac = 1.;
    } else {
	addfac = 2.;
	mulfac = 6.;
    }
    em = (doublereal) (*m);
    en = (doublereal) (*n);
    ek = (doublereal) (*kl);

/*     ---------------------------------   
       GE:  GEneral rectangular matrices   
       --------------------------------- */

    if (lsamen_(&c__2, c2, "GE")) {

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

	if (lsamen_(&c__3, c3, "TRF")) {
	    emn = (doublereal) min(*m,*n);
	    adds = emn * (em * en - (em + en) * (emn + 1.) / 2. + (emn + 1.) *
		     (emn * 2. + 1.) / 6.);
	    mults = adds + emn * (em - (emn + 1.) / 2.);

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

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

/*        xGETRI:  N  =>  M */

	} else if (lsamen_(&c__3, c3, "TRI")) {
	    mults = em * (em * (em * .66666666666666663 + .5) + 
		    .83333333333333337);
	    adds = em * (em * (em * .66666666666666663 - 1.5) + 
		    .83333333333333337);

/*        xGEQRF or xGEQLF:  M, N  =>  M, N */

	} else if (lsamen_(&c__3, c3, "QRF") || lsamen_(
		&c__3, c3, "QR2") || lsamen_(&c__3, c3, 
		"QLF") || lsamen_(&c__3, c3, "QL2")) {
	    if (*m >= *n) {
		mults = en * (em + 3.8333333333333335 + en / 2. + en * (em - 
			en / 3.));
		adds = en * (en * (em - en / 3. + .5) + .83333333333333337);
	    } else {
		mults = em * (en * 2. + 3.8333333333333335 - em / 2. + em * (
			en - em / 3.));
		adds = em * (en + .83333333333333337 - em / 2. + em * (en - 
			em / 3.));
	    }

/*        xGERQF or xGELQF:  M, N  =>  M, N */

	} else if (lsamen_(&c__3, c3, "RQF") || lsamen_(
		&c__3, c3, "RQ2") || lsamen_(&c__3, c3, 
		"LQF") || lsamen_(&c__3, c3, "LQ2")) {
	    if (*m >= *n) {
		mults = en * (em + 4.833333333333333 + en / 2. + en * (em - 
			en / 3.));
		adds = en * (em + .83333333333333337 + en * (em - en / 3. - 
			.5));
	    } else {
		mults = em * (en * 2. + 4.833333333333333 - em / 2. + em * (
			en - em / 3.));
		adds = em * (em / 2. + .83333333333333337 + em * (en - em / 
			3.));
	    }

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

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

/*        xGEQRS or xGERQS:  M, N, NRHS  =>  M, N, KL */

	} else if (lsamen_(&c__3, c3, "QRS") || lsamen_(
		&c__3, c3, "RQS")) {
	    mults = ek * (en * (2. - ek) + em * (en * 2. + (em + 1.) / 2.));
	    adds = ek * (en * (1. - ek) + em * (en * 2. + (em - 1.) / 2.));

/*        xGELQS or xGEQLS:  M, N, NRHS  =>  M, N, KL */

	} else if (lsamen_(&c__3, c3, "LQS") || lsamen_(
		&c__3, c3, "QLS")) {
	    mults = ek * (em * (2. - ek) + en * (em * 2. + (en + 1.) / 2.));
	    adds = ek * (em * (1. - ek) + en * (em * 2. + (en - 1.) / 2.));

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

	} else if (lsamen_(&c__3, c3, "BRD")) {
	    if (*m >= *n) {
		mults = en * (en * (em * 2. - en * .66666666666666663 + 2.) + 
			6.666666666666667);
		adds = en * (en - em + 1.6666666666666667 + en * (em * 2. - 
			en * .66666666666666663));
	    } else {
		mults = em * (em * (en * 2. - em * .66666666666666663 + 2.) + 
			6.666666666666667);
		adds = em * (em - en + 1.6666666666666667 + em * (en * 2. - 
			em * .66666666666666663));
	    }

/*        xGEHRD:  N  =>  M */

	} else if (lsamen_(&c__3, c3, "HRD")) {
	    if (*m == 1) {
		mults = 0.;
		adds = 0.;
	    } else {
		mults = em * (em * (em * 1.6666666666666667 + .5) - 
			1.1666666666666667) - 13.;
		adds = em * (em * (em * 1.6666666666666667 - 1.) - 
			.66666666666666663) - 8.;
	    }

	}

/*     ----------------------------   
       GB:  General Banded matrices   
       ----------------------------   
          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")) {

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

	if (lsamen_(&c__3, c3, "TRF")) {
	    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