/* minfit.f -- translated by f2c (version 19961017).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include <math.h>
#include "f2c.h"

/* Table of constant values */

static doublereal c_b39 = 1.;

/* Subroutine */ int minfit_(integer *nm, integer *m, integer *n, doublereal *
	a, doublereal *w, integer *ip, doublereal *b, integer *ierr, 
	doublereal *rv1)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
    doublereal d__1, d__2, d__3, d__4;

    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    doublereal c__, f, g, h__;
    integer i__, j, k, l=0;
    doublereal s, x, y, z__, scale;
    integer i1, k1, l1=0, m1, ii, kk, ll;
    extern doublereal pythag_(doublereal *, doublereal *);
    integer its;
    doublereal tst1, tst2;



/*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT, */
/*     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH. */
/*     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971). */

/*     THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR */
/*                                                        T */
/*     SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV  OF A REAL */
/*                                         T */
/*     M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U.  HOUSEHOLDER 
*/
/*     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED. */

/*     ON INPUT */

/*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
/*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
/*          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST */
/*          AS LARGE AS THE MAXIMUM OF M AND N. */

/*        M IS THE NUMBER OF ROWS OF A AND B. */

/*        N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V. */

/*        A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM. */

/*        IP IS THE NUMBER OF COLUMNS OF B.  IP CAN BE ZERO. */

/*        B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM */
/*          IF IP IS NOT ZERO.  OTHERWISE B IS NOT REFERENCED. */

/*     ON OUTPUT */

/*        A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE */
/*          DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS.  IF AN */
/*          ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO */
/*          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT. */

/*        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE */
/*          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN */
/*          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT */
/*          FOR INDICES IERR+1,IERR+2,...,N. */

/*                                   T */
/*        B HAS BEEN OVERWRITTEN BY U B.  IF AN ERROR EXIT IS MADE, */
/*                       T */
/*          THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT */
/*          SINGULAR VALUES SHOULD BE CORRECT. */

/*        IERR IS SET TO */
/*          ZERO       FOR NORMAL RETURN, */
/*          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN */
/*                     DETERMINED AFTER 30 ITERATIONS. */

/*        RV1 IS A TEMPORARY STORAGE ARRAY. */

/*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */

/*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
/*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
*/

/*     THIS VERSION DATED AUGUST 1983. */

/*     ------------------------------------------------------------------ 
*/

    /* Parameter adjustments */
    --rv1;
    --w;
    a_dim1 = *nm;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    b_dim1 = *nm;
    b_offset = b_dim1 + 1;
    b -= b_offset;

    /* Function Body */
    *ierr = 0;
/*     .......... HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM .......... */
    g = 0.;
    scale = 0.;
    x = 0.;

    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	l = i__ + 1;
	rv1[i__] = scale * g;
	g = 0.;
	s = 0.;
	scale = 0.;
	if (i__ > *m) {
	    goto L210;
	}

	i__2 = *m;
	for (k = i__; k <= i__2; ++k) {
/* L120: */
	    scale += (d__1 = a[k + i__ * a_dim1], abs(d__1));
	}

	if (scale == 0.) {
	    goto L210;
	}

	i__2 = *m;
	for (k = i__; k <= i__2; ++k) {
	    a[k + i__ * a_dim1] /= scale;
/* Computing 2nd power */
	    d__1 = a[k + i__ * a_dim1];
	    s += d__1 * d__1;
/* L130: */
	}

	f = a[i__ + i__ * a_dim1];
	d__1 = sqrt(s);
	g = -d_sign(&d__1, &f);
	h__ = f * g - s;
	a[i__ + i__ * a_dim1] = f - g;
	if (i__ == *n) {
	    goto L160;
	}

	i__2 = *n;
	for (j = l; j <= i__2; ++j) {
	    s = 0.;

	    i__3 = *m;
	    for (k = i__; k <= i__3; ++k) {
/* L140: */
		s += a[k + i__ * a_dim1] * a[k + j * a_dim1];
	    }

	    f = s / h__;

	    i__3 = *m;
	    for (k = i__; k <= i__3; ++k) {
		a[k + j * a_dim1] += f * a[k + i__ * a_dim1];
/* L150: */
	    }
	}

L160:
	if (*ip == 0) {
	    goto L190;
	}

	i__3 = *ip;
	for (j = 1; j <= i__3; ++j) {
	    s = 0.;

	    i__2 = *m;
	    for (k = i__; k <= i__2; ++k) {
/* L170: */
		s += a[k + i__ * a_dim1] * b[k + j * b_dim1];
	    }

	    f = s / h__;

	    i__2 = *m;
	    for (k = i__; k <= i__2; ++k) {
		b[k + j * b_dim1] += f * a[k + i__ * a_dim1];
/* L180: */
	    }
	}

L190:
	i__2 = *m;
	for (k = i__; k <= i__2; ++k) {
/* L200: */
	    a[k + i__ * a_dim1] = scale * a[k + i__ * a_dim1];
	}

L210:
	w[i__] = scale * g;
	g = 0.;
	s = 0.;
	scale = 0.;
	if (i__ > *m || i__ == *n) {
	    goto L290;
	}

	i__2 = *n;
	for (k = l; k <= i__2; ++k) {
/* L220: */
	    scale += (d__1 = a[i__ + k * a_dim1], abs(d__1));
	}

	if (scale == 0.) {
	    goto L290;
	}

	i__2 = *n;
	for (k = l; k <= i__2; ++k) {
	    a[i__ + k * a_dim1] /= scale;
/* Computing 2nd power */
	    d__1 = a[i__ + k * a_dim1];
	    s += d__1 * d__1;
/* L230: */
	}

	f = a[i__ + l * a_dim1];
	d__1 = sqrt(s);
	g = -d_sign(&d__1, &f);
	h__ = f * g - s;
	a[i__ + l * a_dim1] = f - g;

	i__2 = *n;
	for (k = l; k <= i__2; ++k) {
/* L240: */
	    rv1[k] = a[i__ + k * a_dim1] / h__;
	}

	if (i__ == *m) {
	    goto L270;
	}

	i__2 = *m;
	for (j = l; j <= i__2; ++j) {
	    s = 0.;

	    i__3 = *n;
	    for (k = l; k <= i__3; ++k) {
/* L250: */
		s += a[j + k * a_dim1] * a[i__ + k * a_dim1];
	    }

	    i__3 = *n;
	    for (k = l; k <= i__3; ++k) {
		a[j + k * a_dim1] += s * rv1[k];
/* L260: */
	    }
	}

L270:
	i__3 = *n;
	for (k = l; k <= i__3; ++k) {
/* L280: */
	    a[i__ + k * a_dim1] = scale * a[i__ + k * a_dim1];
	}

L290:
/* Computing MAX */
	d__3 = x, d__4 = (d__1 = w[i__], abs(d__1)) + (d__2 = rv1[i__], abs(
		d__2));
	x = max(d__3,d__4);
/* L300: */
    }
/*     .......... ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS. */
/*                FOR I=N STEP -1 UNTIL 1 DO -- .......... */
    i__1 = *n;
    for (ii = 1; ii <= i__1; ++ii) {
	i__ = *n + 1 - ii;
	if (i__ == *n) {
	    goto L390;
	}
	if (g == 0.) {
	    goto L360;
	}

	i__3 = *n;
	for (j = l; j <= i__3; ++j) {
/*     .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ......
.... */
/* L320: */
	    a[j + i__ * a_dim1] = a[i__ + j * a_dim1] / a[i__ + l * a_dim1] / 
		    g;
	}

	i__3 = *n;
	for (j = l; j <= i__3; ++j) {
	    s = 0.;

	    i__2 = *n;
	    for (k = l; k <= i__2; ++k) {
/* L340: */
		s += a[i__ + k * a_dim1] * a[k + j * a_dim1];
	    }

	    i__2 = *n;
	    for (k = l; k <= i__2; ++k) {
		a[k + j * a_dim1] += s * a[k + i__ * a_dim1];
/* L350: */
	    }
	}

L360:
	i__2 = *n;
	for (j = l; j <= i__2; ++j) {
	    a[i__ + j * a_dim1] = 0.;
	    a[j + i__ * a_dim1] = 0.;
/* L380: */
	}

L390:
	a[i__ + i__ * a_dim1] = 1.;
	g = rv1[i__];
	l = i__;
/* L400: */
    }

    if (*m >= *n || *ip == 0) {
	goto L510;
    }
    m1 = *m + 1;

    i__1 = *n;
    for (i__ = m1; i__ <= i__1; ++i__) {

	i__2 = *ip;
	for (j = 1; j <= i__2; ++j) {
	    b[i__ + j * b_dim1] = 0.;
/* L500: */
	}
    }
/*     .......... DIAGONALIZATION OF THE BIDIAGONAL FORM .......... */
L510:
    tst1 = x;
/*     .......... FOR K=N STEP -1 UNTIL 1 DO -- .......... */
    i__2 = *n;
    for (kk = 1; kk <= i__2; ++kk) {
	k1 = *n - kk;
	k = k1 + 1;
	its = 0;
/*     .......... TEST FOR SPLITTING. */
/*                FOR L=K STEP -1 UNTIL 1 DO -- .......... */
L520:
	i__1 = k;
	for (ll = 1; ll <= i__1; ++ll) {
	    l1 = k - ll;
	    l = l1 + 1;
	    tst2 = tst1 + (d__1 = rv1[l], abs(d__1));
	    if (tst2 == tst1) {
		goto L565;
	    }
/*     .......... RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT */
/*                THROUGH THE BOTTOM OF THE LOOP .......... */
	    tst2 = tst1 + (d__1 = w[l1], abs(d__1));
	    if (tst2 == tst1) {
		goto L540;
	    }
/* L530: */
	}
/*     .......... CANCELLATION OF RV1(L) IF L GREATER THAN 1 .........
. */
L540:
	c__ = 0.;
	s = 1.;

	i__1 = k;
	for (i__ = l; i__ <= i__1; ++i__) {
	    f = s * rv1[i__];
	    rv1[i__] = c__ * rv1[i__];
	    tst2 = tst1 + abs(f);
	    if (tst2 == tst1) {
		goto L565;
	    }
	    g = w[i__];
	    h__ = pythag_(&f, &g);
	    w[i__] = h__;
	    c__ = g / h__;
	    s = -f / h__;
	    if (*ip == 0) {
		goto L560;
	    }

	    i__3 = *ip;
	    for (j = 1; j <= i__3; ++j) {
		y = b[l1 + j * b_dim1];
		z__ = b[i__ + j * b_dim1];
		b[l1 + j * b_dim1] = y * c__ + z__ * s;
		b[i__ + j * b_dim1] = -y * s + z__ * c__;
/* L550: */
	    }

L560:
	    ;
	}
/*     .......... TEST FOR CONVERGENCE .......... */
L565:
	z__ = w[k];
	if (l == k) {
	    goto L650;
	}
/*     .......... SHIFT FROM BOTTOM 2 BY 2 MINOR .......... */
	if (its == 30) {
	    goto L1000;
	}
	++its;
	x = w[l];
	y = w[k1];
	g = rv1[k1];
	h__ = rv1[k];
	f = ((g + z__) / h__ * ((g - z__) / y) + y / h__ - h__ / y) * .5;
	g = pythag_(&f, &c_b39);
	f = x - z__ / x * z__ + h__ / x * (y / (f + d_sign(&g, &f)) - h__);
/*     .......... NEXT QR TRANSFORMATION .......... */
	c__ = 1.;
	s = 1.;

	i__1 = k1;
	for (i1 = l; i1 <= i__1; ++i1) {
	    i__ = i1 + 1;
	    g = rv1[i__];
	    y = w[i__];
	    h__ = s * g;
	    g = c__ * g;
	    z__ = pythag_(&f, &h__);
	    rv1[i1] = z__;
	    c__ = f / z__;
	    s = h__ / z__;
	    f = x * c__ + g * s;
	    g = -x * s + g * c__;
	    h__ = y * s;
	    y *= c__;

	    i__3 = *n;
	    for (j = 1; j <= i__3; ++j) {
		x = a[j + i1 * a_dim1];
		z__ = a[j + i__ * a_dim1];
		a[j + i1 * a_dim1] = x * c__ + z__ * s;
		a[j + i__ * a_dim1] = -x * s + z__ * c__;
/* L570: */
	    }

	    z__ = pythag_(&f, &h__);
	    w[i1] = z__;
/*     .......... ROTATION CAN BE ARBITRARY IF Z IS ZERO .........
. */
	    if (z__ == 0.) {
		goto L580;
	    }
	    c__ = f / z__;
	    s = h__ / z__;
L580:
	    f = c__ * g + s * y;
	    x = -s * g + c__ * y;
	    if (*ip == 0) {
		goto L600;
	    }

	    i__3 = *ip;
	    for (j = 1; j <= i__3; ++j) {
		y = b[i1 + j * b_dim1];
		z__ = b[i__ + j * b_dim1];
		b[i1 + j * b_dim1] = y * c__ + z__ * s;
		b[i__ + j * b_dim1] = -y * s + z__ * c__;
/* L590: */
	    }

L600:
	    ;
	}

	rv1[l] = 0.;
	rv1[k] = f;
	w[k] = x;
	goto L520;
/*     .......... CONVERGENCE .......... */
L650:
	if (z__ >= 0.) {
	    goto L700;
	}
/*     .......... W(K) IS MADE NON-NEGATIVE .......... */
	w[k] = -z__;

	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
/* L690: */
	    a[j + k * a_dim1] = -a[j + k * a_dim1];
	}

L700:
	;
    }

    goto L1001;
/*     .......... SET ERROR -- NO CONVERGENCE TO A */
/*                SINGULAR VALUE AFTER 30 ITERATIONS .......... */
L1000:
    *ierr = k;
L1001:
    return 0;
} /* minfit_ */

