sggsvp.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*     Determine the effective rank of matrix B. */

/*<       L = 0 >*/
    *l = 0;
/*<       DO 20 I = 1, MIN( P, N ) >*/
    i__1 = min(*p,*n);
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<    >*/
        if ((r__1 = b[i__ + i__ * b_dim1], dabs(r__1)) > *tolb) {
            ++(*l);
        }
/*<    20 CONTINUE >*/
/* L20: */
    }

/*<       IF( WANTV ) THEN >*/
    if (wantv) {

/*        Copy the details of V, and form V. */

/*<          CALL SLASET( 'Full', P, P, ZERO, ZERO, V, LDV ) >*/
        slaset_("Full", p, p, &c_b12, &c_b12, &v[v_offset], ldv, (ftnlen)4);
/*<    >*/
        if (*p > 1) {
            i__1 = *p - 1;
            slacpy_("Lower", &i__1, n, &b[b_dim1 + 2], ldb, &v[v_dim1 + 2], 
                    ldv, (ftnlen)5);
        }
/*<          CALL SORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) >*/
        i__1 = min(*p,*n);
        sorg2r_(p, p, &i__1, &v[v_offset], ldv, &tau[1], &work[1], info);
/*<       END IF >*/
    }

/*     Clean up B */

/*<       DO 40 J = 1, L - 1 >*/
    i__1 = *l - 1;
    for (j = 1; j <= i__1; ++j) {
/*<          DO 30 I = J + 1, L >*/
        i__2 = *l;
        for (i__ = j + 1; i__ <= i__2; ++i__) {
/*<             B( I, J ) = ZERO >*/
            b[i__ + j * b_dim1] = (float)0.;
/*<    30    CONTINUE >*/
/* L30: */
        }
/*<    40 CONTINUE >*/
/* L40: */
    }
/*<    >*/
    if (*p > *l) {
        i__1 = *p - *l;
        slaset_("Full", &i__1, n, &c_b12, &c_b12, &b[*l + 1 + b_dim1], ldb, (
                ftnlen)4);
    }

/*<       IF( WANTQ ) THEN >*/
    if (wantq) {

/*        Set Q = I and Update Q := Q*P */

/*<          CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) >*/
        slaset_("Full", n, n, &c_b12, &c_b22, &q[q_offset], ldq, (ftnlen)4);
/*<          CALL SLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) >*/
        slapmt_(&forwrd, n, n, &q[q_offset], ldq, &iwork[1]);
/*<       END IF >*/
    }

/*<       IF( P.GE.L .AND. N.NE.L ) THEN >*/
    if (*p >= *l && *n != *l) {

/*        RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z */

/*<          CALL SGERQ2( L, N, B, LDB, TAU, WORK, INFO ) >*/
        sgerq2_(l, n, &b[b_offset], ldb, &tau[1], &work[1], info);

/*        Update A := A*Z' */

/*<    >*/
        sormr2_("Right", "Transpose", m, n, l, &b[b_offset], ldb, &tau[1], &a[
                a_offset], lda, &work[1], info, (ftnlen)5, (ftnlen)9);

/*<          IF( WANTQ ) THEN >*/
        if (wantq) {

/*           Update Q := Q*Z' */

/*<    >*/
            sormr2_("Right", "Transpose", n, n, l, &b[b_offset], ldb, &tau[1],
                     &q[q_offset], ldq, &work[1], info, (ftnlen)5, (ftnlen)9);
/*<          END IF >*/
        }

/*        Clean up B */

/*<          CALL SLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB ) >*/
        i__1 = *n - *l;
        slaset_("Full", l, &i__1, &c_b12, &c_b12, &b[b_offset], ldb, (ftnlen)
                4);
/*<          DO 60 J = N - L + 1, N >*/
        i__1 = *n;
        for (j = *n - *l + 1; j <= i__1; ++j) {
/*<             DO 50 I = J - N + L + 1, L >*/
            i__2 = *l;
            for (i__ = j - *n + *l + 1; i__ <= i__2; ++i__) {
/*<                B( I, J ) = ZERO >*/
                b[i__ + j * b_dim1] = (float)0.;
/*<    50       CONTINUE >*/
/* L50: */
            }
/*<    60    CONTINUE >*/
/* L60: */
        }

/*<       END IF >*/
    }

/*     Let              N-L     L */
/*                A = ( A11    A12 ) M, */

/*     then the following does the complete QR decomposition of A11: */

/*              A11 = U*(  0  T12 )*P1' */
/*                      (  0   0  ) */

/*<       DO 70 I = 1, N - L >*/
    i__1 = *n - *l;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          IWORK( I ) = 0 >*/
        iwork[i__] = 0;
/*<    70 CONTINUE >*/
/* L70: */
    }
/*<       CALL SGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, INFO ) >*/
    i__1 = *n - *l;
    sgeqpf_(m, &i__1, &a[a_offset], lda, &iwork[1], &tau[1], &work[1], info);

/*     Determine the effective rank of A11 */

/*<       K = 0 >*/
    *k = 0;
/*<       DO 80 I = 1, MIN( M, N-L ) >*/
/* Computing MIN */
    i__2 = *m, i__3 = *n - *l;
    i__1 = min(i__2,i__3);
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<    >*/
        if ((r__1 = a[i__ + i__ * a_dim1], dabs(r__1)) > *tola) {
            ++(*k);
        }
/*<    80 CONTINUE >*/
/* L80: */
    }

/*     Update A12 := U'*A12, where A12 = A( 1:M, N-L+1:N ) */

/*<    >*/
/* Computing MIN */
    i__2 = *m, i__3 = *n - *l;
    i__1 = min(i__2,i__3);
    sorm2r_("Left", "Transpose", m, l, &i__1, &a[a_offset], lda, &tau[1], &a[(
            *n - *l + 1) * a_dim1 + 1], lda, &work[1], info, (ftnlen)4, (
            ftnlen)9);

/*<       IF( WANTU ) THEN >*/
    if (wantu) {

/*        Copy the details of U, and form U */

/*<          CALL SLASET( 'Full', M, M, ZERO, ZERO, U, LDU ) >*/
        slaset_("Full", m, m, &c_b12, &c_b12, &u[u_offset], ldu, (ftnlen)4);
/*<    >*/
        if (*m > 1) {
            i__1 = *m - 1;
            i__2 = *n - *l;
            slacpy_("Lower", &i__1, &i__2, &a[a_dim1 + 2], lda, &u[u_dim1 + 2]
                    , ldu, (ftnlen)5);
        }
/*<          CALL SORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) >*/
/* Computing MIN */
        i__2 = *m, i__3 = *n - *l;
        i__1 = min(i__2,i__3);
        sorg2r_(m, m, &i__1, &u[u_offset], ldu, &tau[1], &work[1], info);
/*<       END IF >*/
    }

/*<       IF( WANTQ ) THEN >*/
    if (wantq) {

/*        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1 */

/*<          CALL SLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) >*/
        i__1 = *n - *l;
        slapmt_(&forwrd, n, &i__1, &q[q_offset], ldq, &iwork[1]);
/*<       END IF >*/
    }

/*     Clean up A: set the strictly lower triangular part of */
/*     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. */

/*<       DO 100 J = 1, K - 1 >*/
    i__1 = *k - 1;
    for (j = 1; j <= i__1; ++j) {
/*<          DO 90 I = J + 1, K >*/
        i__2 = *k;
        for (i__ = j + 1; i__ <= i__2; ++i__) {
/*<             A( I, J ) = ZERO >*/
            a[i__ + j * a_dim1] = (float)0.;
/*<    90    CONTINUE >*/
/* L90: */
        }
/*<   100 CONTINUE >*/
/* L100: */
    }
/*<    >*/
    if (*m > *k) {
        i__1 = *m - *k;
        i__2 = *n - *l;
        slaset_("Full", &i__1, &i__2, &c_b12, &c_b12, &a[*k + 1 + a_dim1], 
                lda, (ftnlen)4);
    }

/*<       IF( N-L.GT.K ) THEN >*/
    if (*n - *l > *k) {

/*        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 */

/*<          CALL SGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) >*/
        i__1 = *n - *l;
        sgerq2_(k, &i__1, &a[a_offset], lda, &tau[1], &work[1], info);

/*<          IF( WANTQ ) THEN >*/
        if (wantq) {

/*           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1' */

/*<    >*/
            i__1 = *n - *l;
            sormr2_("Right", "Transpose", n, &i__1, k, &a[a_offset], lda, &
                    tau[1], &q[q_offset], ldq, &work[1], info, (ftnlen)5, (
                    ftnlen)9);
/*<          END IF >*/
        }

/*        Clean up A */

/*<          CALL SLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA ) >*/
        i__1 = *n - *l - *k;
        slaset_("Full", k, &i__1, &c_b12, &c_b12, &a[a_offset], lda, (ftnlen)
                4);
/*<          DO 120 J = N - L - K + 1, N - L >*/
        i__1 = *n - *l;
        for (j = *n - *l - *k + 1; j <= i__1; ++j) {
/*<             DO 110 I = J - N + L + K + 1, K >*/
            i__2 = *k;
            for (i__ = j - *n + *l + *k + 1; i__ <= i__2; ++i__) {
/*<                A( I, J ) = ZERO >*/
                a[i__ + j * a_dim1] = (float)0.;
/*<   110       CONTINUE >*/
/* L110: */
            }
/*<   120    CONTINUE >*/
/* L120: */
        }

/*<       END IF >*/
    }

/*<       IF( M.GT.K ) THEN >*/
    if (*m > *k) {

/*        QR factorization of A( K+1:M,N-L+1:N ) */

/*<          CALL SGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) >*/
        i__1 = *m - *k;
        sgeqr2_(&i__1, l, &a[*k + 1 + (*n - *l + 1) * a_dim1], lda, &tau[1], &
                work[1], info);

/*<          IF( WANTU ) THEN >*/
        if (wantu) {

/*           Update U(:,K+1:M) := U(:,K+1:M)*U1 */

/*<    >*/
            i__1 = *m - *k;
/* Computing MIN */
            i__3 = *m - *k;
            i__2 = min(i__3,*l);
            sorm2r_("Right", "No transpose", m, &i__1, &i__2, &a[*k + 1 + (*n 
                    - *l + 1) * a_dim1], lda, &tau[1], &u[(*k + 1) * u_dim1 + 
                    1], ldu, &work[1], info, (ftnlen)5, (ftnlen)12);
/*<          END IF >*/
        }

/*        Clean up */

/*<          DO 140 J = N - L + 1, N >*/
        i__1 = *n;
        for (j = *n - *l + 1; j <= i__1; ++j) {
/*<             DO 130 I = J - N + K + L + 1, M >*/
            i__2 = *m;
            for (i__ = j - *n + *k + *l + 1; i__ <= i__2; ++i__) {
/*<                A( I, J ) = ZERO >*/
                a[i__ + j * a_dim1] = (float)0.;
/*<   130       CONTINUE >*/
/* L130: */
            }
/*<   140    CONTINUE >*/
/* L140: */
        }

/*<       END IF >*/
    }

/*<       RETURN >*/
    return 0;

/*     End of SGGSVP */

/*<       END >*/
} /* sggsvp_ */

#ifdef __cplusplus
        }
#endif