/*

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

Authors

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

Per-Gunnar Martinsson

  Dept. of Applied Mathematics,

  University of Colorado at Boulder,

  526 UCB, Boulder, CO 80309-0526, USA

Gregorio Quintana-Orti

  Depto. de Ingenieria y Ciencia de Computadores,

  Universitat Jaume I,

  12.071 Castellon, Spain

Nathan Heavner

  Dept. of Applied Mathematics,

  University of Colorado at Boulder,

  526 UCB, Boulder, CO 80309-0526, USA

Robert van de Geijn

  Dept. of Computer Science and Institute for Computational Engineering and

  Sciences,

  The University of Texas at Austin

  Austin, TX.

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

Copyright

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

Copyright (C) 2016,

  Universitat Jaume I,

  University of Colorado at Boulder,

  The University of Texas at Austin.

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

Disclaimer

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

This code is distributed in the hope that it will be useful, but

WITHOUT ANY WARRANTY EXPRESSED OR IMPLIED.

*/

#include <omp.h>

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

#include "blas_lapack_prototypes.h"

#include "dgeqp4.h"

// Matrices with dimensions smaller than THRESHOLD_FOR_DGEQPF are processed

// with LAPACK's routine dgeqpf.

// Matrices with dimensions between THRESHOLD_FOR_DGEQPF and

// THRESHOLD_FOR_DGEQP3 are processed with LAPACK's routine dgeqp3.

// Matrices with dimensions larger than THRESHOLD_FOR_DGEQP3 are processed

// with the new HQRRP code.

#define THRESHOLD_FOR_DGEQPF    0

#define THRESHOLD_FOR_DGEQP3  500

// ============================================================================

// Definition of macros.

#define max( a, b )  ( (a) > (b) ? (a) : (b) )

#define min( a, b )  ( (a) > (b) ? (b) : (a) )

#define dabs( a )    ( (a) >= 0.0 ? (a) : -(a) )

// ============================================================================

// Compilation declarations.

#undef CHECK_DOWNDATING_OF_Y

// ============================================================================

// Declaration of local prototypes.

static int dgeqp4_Normal_random_matrix( int m_A, int n_A,

               double * buff_A, int ldim_A );

static double dgeqp4_Normal_random_number( double mu, double sigma );

static int dgeqp4_Downdate_Y(

               int m_U11, int n_U11, double * buff_U11, int ldim_U11,

               int m_U21, int n_U21, double * buff_U21, int ldim_U21,

               int m_A12, int n_A12, double * buff_A12, int ldim_A12,

               int m_T, int n_T, double * buff_T, int ldim_T,

               int m_Y2, int n_Y2, double * buff_Y2, int ldim_Y2,

               int m_G1, int n_G1, double * buff_G1, int ldim_G1,

               int m_G2, int n_G2, double * buff_G2, int ldim_G2 );

static int dgeqp4_Apply_Q_WY_lhfc_blk_var4(

               int m_U, int n_U, double * buff_U, int ldim_U,

               int m_T, int n_T, double * buff_T, int ldim_T,

               int m_B, int n_B, double * buff_B, int ldim_B );

static int dgeqp4_Apply_Q_WY_rnfc_blk_var4(

               int m_U, int n_U, double * buff_U, int ldim_U,

               int m_T, int n_T, double * buff_T, int ldim_T,

               int m_B, int n_B, double * buff_B, int ldim_B );

static int dgeqp4_QRPmod_WY_unb_var4( int pivoting, int num_stages,

               int m_A, int n_A, double * buff_A, int ldim_A,

               int * buff_p, double * buff_t,

               int pivot_B, int m_B, double * buff_B, int ldim_B,

               int pivot_C, int m_C, double * buff_C, int ldim_C,

               int build_T, double * buff_T, int ldim_T );

static int dgeqp4_QRP_compute_norms(

               int m_A, int n_A, double * buff_A, int ldim_A,

               double * buff_d, double * buff_e );

static int dgeqp4_QRP_downdate_partial_norms( int m_A, int n_A,

               double * buff_d,  int st_d,

               double * buff_e,  int st_e,

               double * buff_wt, int st_wt,

               double * buff_A,  int ldim_A );

static int dgeqp4_QRP_pivot_G_B_C( int j_max_col,

               int m_G, double * buff_G, int ldim_G,

               int pivot_B, int m_B, double * buff_B, int ldim_B,

               int pivot_C, int m_C, double * buff_C, int ldim_C,

               int * buff_p,

               double * buff_d, double * buff_e );

// ============================================================================

        double * work, int * lwork, int * info ) {

//

// This routine is plug compatible with LAPACK's routine dgeqp3.

// It computes the new HQRRP while keeping the same header as LAPACK's dgeqp3.

// It uses dgeqpf or dgeqp3 for small matrices. The thresholds are defined in

// constants THRESHOLD_FOR_DGEQPF and THRESHOLD_FOR_DGEQP3.

//

          m_A, n_A, mn_A, ldim_A, lquery, nb, num_factorized_fixed_cols,

          minus_info, iws, lwkopt, j, k, num_fixed_cols, n_rest, itmp;

  int     * previous_jpvt;

  int     ilaenv_();

  // Some initializations.

  // Check input arguments.

  }

      iws    = 1;

      lwkopt = 1;

    } else {

                        & i_minus_one );

    }

    }

  }

    return;

  }

  // Quick return if possible.

    return;

  }

  // Use LAPACK's DGEQPF or DGEQP3 for small matrices.

  //if( mn_A < THRESHOLD_FOR_DGEQPF ) {

  //  // Call to LAPACK routine.

  //  //// printf( "Calling dgeqpf\n" );

  //  dgeqpf_( m, n, A, lda, jpvt, tau, work, info );

  //  return;

  //} else

  // Use LAPACK's DGEQP3 for small matrices.

    //// printf( "Calling dgeqp3\n" );

  }

  // Move initial columns up front.

  num_fixed_cols = 0;

        //// printf( "Swapping columns: %d %d \n", j, num_fixed_cols );

      } else {

      }

    } else {

    }

  }

  // Factorize fixed columns at the front.

             info );

    }

               & m_A, & n_rest, & num_factorized_fixed_cols,

               A, & ldim_A, tau,

               work, lwork, info );

      }

    }

  }

  // Create intermediate jpvt vector.

  // Save a copy of jpvt vector.

    // Copy vector.

    }

  }

  // Factorize free columns at the bottom with default values:

  // nb_alg = 64, pp = 10, panel_pivoting = 1.

        m_A - num_factorized_fixed_cols, n_A - num_factorized_fixed_cols,

            ldim_A,

        64, 10, 1 );

  }

  // Pivot block above factorized block by dgeqp4_HQRRP.

    // Pivot block above factorized block.

      //// printf( "%% Processing j: %d \n", j );

          //// printf( "%%   Found j: %d  k: %d \n", j, k );

          break;

        }

      }

      // Swap vector previous_jpvt and block above factorized block.

        // Swap elements in previous_jpvt.

        //// printf( "%%   Swapping  j: %d  k: %d \n", j, k );

        // Swap columns in block above factorized block.

      }

    }

  }

  // Remove intermediate jpvt vector.

  // Return workspace length required.

}

// ============================================================================

        int * buff_jpvt, double * buff_tau,

        int nb_alg, int pp, int panel_pivoting ) {

//

// HQRRP: It computes the Householder QR with Randomized Pivoting of matrix A.

// This routine is almost compatible with LAPACK's dgeqp3.

// The main difference is that this routine does not manage fixed columns.

//

// Main features:

//   * BLAS-3 based.

//   * Norm downdating method by Drmac.

//   * Downdating for computing Y.

//   * No use of libflame.

//   * Compact WY transformations are used instead of UT transformations.

//   * LAPACK's routine dlarfb is used to apply block transformations.

//

// Arguments:

// ----------

// m_A:            Number of rows of matrix A.

// n_A:            Number of columns of matrix A.

// buff_A:         Address/pointer of/to data in matrix A. Matrix A must be

//                 stored in column-order.

// ldim_A:         Leading dimension of matrix A.

// buff_jpvt:      Input/output vector with the pivots.

// buff_tau:       Output vector with the tau values of the Householder factors.

// nb_alg:         Block size.

//                 Usual values for nb_alg are 32, 64, etc.

// pp:             Oversampling size.

//                 Usual values for pp are 5, 10, etc.

// panel_pivoting: If panel_pivoting==1, QR with pivoting is applied to

//                 factorize the panels of matrix A. Otherwise, QR without

//                 pivoting is used. Usual value for panel_pivoting is 1.

// Final comments:

// ---------------

// This code has been created from a libflame code. Hence, you can find some

// commented calls to libflame routines. We have left them to make it easier

// to interpret the meaning of the C code.

//

  int     b, j, last_iter, mn_A, m_Y, n_Y, ldim_Y, m_V, n_V, ldim_V,

          m_W, n_W, ldim_W, n_VR, m_AB1, n_AB1, ldim_T1_T,

          m_A11, n_A11, m_A12, n_A12, m_A21, n_A21, m_A22,

          m_G, n_G, ldim_G;

  double  * buff_Y, * buff_V, * buff_W, * buff_VR, * buff_YR,

          * buff_s, * buff_sB, * buff_s1,

          * buff_AR, * buff_AB1, * buff_A01, * buff_Y1, * buff_T1_T,

          * buff_A11, * buff_A21, * buff_A12,

          * buff_Y2, * buff_G, * buff_G1, * buff_G2;

  int     * buff_p, * buff_pB, * buff_p1;

  double down_t = 0.0, unb_t = 0.0, lhfc_t = 0.0;

  // Executable Statements.

  //// printf( "%% dgeqp4_HQRRP_WY_blk_var4.\n" );

  // Check arguments.

             "ERROR in dgeqp4_HQRRP_WY_blk_var4: m_A is < 0.\n" );

             "ERROR in dgeqp4_HQRRP_WY_blk_var4: n_A is < 0.\n" );

             "ERROR in dgeqp4_HQRRP_WY_blk_var4: ldim_A is < max( 1, m_A ).\n" );

  }

  // Some initializations.

  buff_p = buff_jpvt;

  buff_s = buff_tau;

  // Quick return.

    return 0;

  }

  // Initialize the seed for the generator of random numbers.

  // Create auxiliary objects.

  n_V     = n_A;

  m_W     = nb_alg;

  n_W     = n_A;

  ldim_W  = m_W;

  m_G     = nb_alg + pp;

  // Initialize matrices G and Y.

  //// FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,

  ////           FLA_ONE, G, A, FLA_ZERO, Y );

          & d_one, buff_G, & ldim_G, buff_A, & ldim_A,

          & d_zero, buff_Y, & ldim_Y );

  //printf( "rand_t %lf gemm_proj_t %lf\n", rand_t, gemm_proj_t ); fflush( stdout );

  // Main Loop.

    // Check whether it is the last iteration.

    // Some initializations for the iteration of this loop.

    m_AB1     = m_A - j;

    n_AB1     = b;

    buff_p1   = & buff_p[ j ];

    buff_s1   = & buff_s[ j ];

    buff_A01  = & buff_A[ 0 + j * ldim_A ];

    buff_Y1   = & buff_Y[ 0 + j * ldim_Y ];

    ldim_T1_T = ldim_W;

    buff_A11 = & buff_A[ j + j * ldim_A ];

    m_A11 = b;

    n_A11 = b;

    n_A21 = b;

    m_A12 = b;

    //// buff_A22 = & buff_A[ min( m_A - 1, j + b ) +

    ////                      min( n_A - 1, j + b ) * ldim_A ];

    m_A22 = max( 0, m_A - j - b );

    //// n_A22 = max( 0, n_A - j - b );

#ifdef CHECK_DOWNDATING_OF_Y

    // Check downdating of matrix Y: Compare downdated matrix Y with

    // matrix Y computed from scratch.

    int     m_cyr, n_cyr, ldim_cyr, m_ABR, ii, jj;

    double  * buff_cyr, aux, sum;

    m_cyr    = m_Y;

    n_cyr    = n_Y - j;

    ldim_cyr = m_cyr;

    m_ABR    = m_A - j;

    buff_cyr = ( double * ) malloc( m_cyr * n_cyr * sizeof( double ) );

    //// FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,

    ////           FLA_ONE, GR, ABR, FLA_ZERO, CYR );

    dgemm_( "No tranpose", "No transpose", & m_cyr, & n_cyr, & m_ABR,

            & d_one, & buff_G[ 0 + j * ldim_G ], & ldim_G,

                     & buff_A[ j + j * ldim_A ], & ldim_A,

            & d_zero, & buff_cyr[ 0 + 0 * ldim_cyr ], & ldim_cyr );

    //// print_double_matrix( "cyr", m_cyr, n_cyr, buff_cyr, ldim_cyr );

    //// print_double_matrix( "y", m_Y, n_Y, buff_Y, ldim_Y );

    sum = 0.0;

    for( jj = 0; jj < n_cyr; jj++ ) {

      for( ii = 0; ii < m_cyr; ii++ ) {

        aux = buff_Y[ ii + ( j + jj ) * ldim_Y ] -

              buff_cyr[ ii + jj * ldim_cyr ];

        sum += aux * aux;

      }

    }

    sum = sqrt( sum );

    printf( "%%  diff between Y and downdated Y: %le\n", sum );

    free( buff_cyr );

#endif

      // Compute QRP of YR, and apply permutations to matrix AR.

      // A copy of YR is made into VR, and permutations are applied to YR.

      //// FLA_Merge_2x1( ATR,

      ////                ABR,   & AR );

      //// FLA_Copy( YR, VR );

      //// FLA_QRPmod_WY_unb_var4( 1, bRow, VR, pB, sB, 1, AR, 1, YR, 0, None );

                                     buff_VR, & ldim_V );

          m_V, n_VR, buff_VR, ldim_V, buff_pB, buff_sB,

          1, m_A, buff_AR, ldim_A,

          1, m_Y, buff_YR, ldim_Y,

          0, buff_Y, ldim_Y );

    }

    //

    // Compute QRP of panel AB1 = [ A11; A21 ].

    // Apply same permutations to A01 and Y1, and build T1_T.

    //

    //// FLA_Part_2x1( W1,   & T1_T,

    ////                     & None,    b, FLA_TOP );

    //// FLA_Merge_2x1( A11,

    ////                A21,   & AB1 );

    //// FLA_QRPmod_WY_unb_var4( panel_pivoting, -1, AB1, p1, s1,

    ////                         1, A01, 1, Y1, 1, T1_T );

        m_AB1, n_AB1, buff_AB1, ldim_A, buff_p1, buff_s1,

        1, j, buff_A01, ldim_A,

        1, m_Y, buff_Y1, ldim_Y,

        1, buff_T1_T, ldim_W );

    //

    // Update the rest of the matrix.

    //

      // Apply the Householder transforms associated with AB1 = [ A11; A21 ]

      // and T1_T to [ A12; A22 ]:

      //   / A12 \ := QB1' / A12 \

      //   \ A22 /         \ A22 /

      // where QB1 is formed from AB1 and T1_T.

      //// MyFLA_Apply_Q_WY_lhfc_blk_var4( A11, A21, T1_T, A12, A22 );

          m_A11 + m_A21, n_A11, buff_A11, ldim_A,

          b, b, buff_T1_T, ldim_W,

          m_A12 + m_A22, n_A12, buff_A12, ldim_A );

    }

    //

    // Downdate matrix Y.

    //

      //// MyFLA_Downdate_Y( A11, A21, A12, T1_T, Y2, G1, G2 );

          m_A11, n_A11, buff_A11, ldim_A,

          m_A21, n_A21, buff_A21, ldim_A,

          m_A12, n_A12, buff_A12, ldim_A,

          b, b, buff_T1_T, ldim_T1_T,

          m_G, b, buff_G1, ldim_G,

    }

  }

  //printf( "down_t %lf, unb_t %lf, lhfc %lf\n", down_t, unb_t, lhfc_t ); fflush( stdout );

  // Remove auxiliary objects.

  //// FLA_Obj_free( & G );

  //// FLA_Obj_free( & Y );

  //// FLA_Obj_free( & V );

  //// FLA_Obj_free( & W );

}

// ============================================================================

               double * buff_A, int ldim_A ) {

//

// It generates a random matrix with normal distribution.

//

  int  i, j;

  // Main loop.

    }

  }

}

// ============================================================================

//

// It computes and returns a normal random number.

//

  static int     alternate_calls = 0;

  static double  b1, b2;

  double         c1, c2, a, factor;

  // Quick return.

  }

  // Main loop.

  do {

}

// ============================================================================

               int m_U11, int n_U11, double * buff_U11, int ldim_U11,

               int m_U21, int n_U21, double * buff_U21, int ldim_U21,

               int m_A12, int n_A12, double * buff_A12, int ldim_A12,

               int m_T, int n_T, double * buff_T, int ldim_T,

               int m_Y2, int n_Y2, double * buff_Y2, int ldim_Y2,

               int m_G1, int n_G1, double * buff_G1, int ldim_G1,

               int m_G2, int n_G2, double * buff_G2, int ldim_G2 ) {

//

// It downdates matrix Y, and updates matrix G.

// Only Y2 of Y is updated.

// Only G1 and G2 of G are updated.

//

// Y2 = Y2 - ( G1 - ( G1*U11 + G2*U21 ) * inv( T11 ) * U11' ) * R12.

//

  int    i, j;

  double * buff_B;

  // Create object B.

  //// FLA_Obj_create_conf_to( FLA_NO_TRANSPOSE, G1, & B );

  // B = G1.

  //// FLA_Copy( G1, B );

                                  buff_B, & ldim_B );

  // B = B * U11.

  //// FLA_Trmm( FLA_RIGHT, FLA_LOWER_TRIANGULAR,

  ////           FLA_NO_TRANSPOSE, FLA_UNIT_DIAG,

  ////           FLA_ONE, U11, B );

          & d_one, buff_U11, & ldim_U11, buff_B, & ldim_B );

  // B = B + G2 * U21.

  //// FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,

  ////           FLA_ONE, G2, U21, FLA_ONE, B );

          & d_one, buff_G2, & ldim_G2, buff_U21, & ldim_U21,

          & d_one, buff_B, & ldim_B );

  // B = B * inv( T11 ).

  //// FLA_Trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,

  ////           FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,

  ////           FLA_ONE, T, B );

  //// dtrsm_( "Right", "Upper", "No transpose", "Non-unit", & m_B, & n_B,

  ////         & d_one, buff_T, & ldim_T, buff_B, & ldim_B );

  // Used dtrmm instead of dtrsm because of using compact WY instead of UT.

          & d_one, buff_T, & ldim_T, buff_B, & ldim_B );

  // B = - B * U11^H.

  //// FLA_Trmm( FLA_RIGHT, FLA_LOWER_TRIANGULAR,

  ////           FLA_CONJ_TRANSPOSE, FLA_UNIT_DIAG,

  ////           FLA_MINUS_ONE, U11, B );

          & d_minus_one, buff_U11, & ldim_U11, buff_B, & ldim_B );

  // B = G1 + B.

  //// FLA_Axpy( FLA_ONE, G1, B );

    }

  }

  // Y2 = Y2 - B * R12.

  //// FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,

  ////           FLA_MINUS_ONE, B, A12, FLA_ONE, Y2 );

          & d_minus_one, buff_B, & ldim_B, buff_A12, & ldim_A12,

          & d_one, buff_Y2, & ldim_Y2 );

  //

  // GR = GR * Q

  //

          m_U11 + m_U21, n_U11, buff_U11, ldim_U11,

          m_T, n_T, buff_T, ldim_T,

          m_G1, n_G1 + n_G2, buff_G1, ldim_G1 );

  // Remove object B.

  //// FLA_Obj_free( & B );

}

// ============================================================================

               int m_U, int n_U, double * buff_U, int ldim_U,

               int m_T, int n_T, double * buff_T, int ldim_T,

               int m_B, int n_B, double * buff_B, int ldim_B ) {

//

// It applies the transpose of a block transformation Q to a matrix B from

// the left:

//   B := Q' * B

// where:

//   Q = I - U * T' * U'.

//

  double  * buff_W;

  int     ldim_W;

  // Create auxiliary object.

  //// FLA_Obj_create_conf_to( FLA_NO_TRANSPOSE, B1, & W );

  // Apply the block transformation.

           & m_B, & n_B, & n_U, buff_U, & ldim_U, buff_T, & ldim_T,

           buff_B, & ldim_B, buff_W, & ldim_W );

  // Remove auxiliary object.

  //// FLA_Obj_free( & W );

}

// ============================================================================

               int m_U, int n_U, double * buff_U, int ldim_U,

               int m_T, int n_T, double * buff_T, int ldim_T,

               int m_B, int n_B, double * buff_B, int ldim_B ) {

//

// It applies a block transformation Q to a matrix B from the right:

//   B = B * Q

// where:

//   Q = I - U * T' * U'.

//

  double  * buff_W;

  int     ldim_W;

  // Create auxiliary object.

  //// FLA_Obj_create_conf_to( FLA_TRANSPOSE, B1, & W );

  // Apply the block transformation.

           & m_B, & n_B, & n_U, buff_U, & ldim_U, buff_T, & ldim_T,

           buff_B, & ldim_B, buff_W, & ldim_W );

  // Remove auxiliary object.

  //// FLA_Obj_free( & W );

}

// ============================================================================

               int m_A, int n_A, double * buff_A, int ldim_A,

               int * buff_p, double * buff_t,

               int pivot_B, int m_B, double * buff_B, int ldim_B,

               int pivot_C, int m_C, double * buff_C, int ldim_C,

               int build_T, double * buff_T, int ldim_T ) {

//

// It computes an unblocked QR factorization of matrix A with or without

// pivoting. Matrices B and C are optionally pivoted, and matrix T is

// optionally built.

//

// Arguments:

// "pivoting": If pivoting==1, then QR factorization with pivoting is used.

// "numstages": It tells the number of columns that are factorized.

//   If "num_stages" is negative, the whole matrix A is factorized.

//   If "num_stages" is positive, only the first "num_stages" are factorized.

// "pivot_B": if "pivot_B" is true, matrix "B" is pivoted too.

// "pivot_C": if "pivot_C" is true, matrix "C" is pivoted too.

// "build_T": if "build_T" is true, matrix "T" is built.

//

  int     j, mn_A, m_a21, m_A22, n_A22, n_dB, idx_max_col,

  double  * buff_d, * buff_e, * buff_workspace, diag;

  int     idamax_();

  //// printf( "dgeqp4_QRPmod_WY_unb_var4. pivoting: %d \n", pivoting );

  // Some initializations.

  // Set the number of stages, if needed.

  }

  // Create auxiliary vectors.

    // Compute initial norms of A into d and e.

  }

  // Main Loop.

    m_A22 = m_A - j - 1;

      // Obtain the index of the column with largest 2-norm.

      // Swap columns of A, B, C, pivots, and norms vectors.

          & buff_d[ j ],

          & buff_e[ j ] );

    }

    // Compute tau1 and u21 from alpha11 and a21 such that tau1 and u21

    // determine a Householder transform H such that applying H from the

    // left to the column vector consisting of alpha11 and a21 annihilates

    // the entries in a21 (and updates alpha11).

    // / a12t \ =  H / a12t \

    // \ A22  /      \ A22  /

    //

    // where H is formed from tau1 and u21.

        & buff_A[ j + j * ldim_A ], & i_one,

        & buff_t[ j ],

        buff_workspace );

      // Update partial column norms.

          & buff_d[ j+1 ], 1,

          & buff_e[ j+1 ], 1,

    }

  }

  // Build T.

             buff_t, buff_T, & ldim_T );

  }

  // Remove auxiliary vectors.

}

// ============================================================================

               int m_A, int n_A, double * buff_A, int ldim_A,

               double * buff_d, double * buff_e ) {

//

// It computes the column norms of matrix A. The norms are stored into

// vectors d and e.

//

  double  dnrm2_();

  // Main loop.

  }

}

// ============================================================================

               double * buff_d,  int st_d,

               double * buff_e,  int st_e,

               double * buff_wt, int st_wt,

               double * buff_A,  int ldim_A ) {

//

// It updates (downdates) the column norms of matrix A. It uses Drmac's method.

//

  double  * ptr_d, * ptr_e, * ptr_wt, * ptr_A;

  double  temp, temp2, temp5, tol3z;

  double  dnrm2_(), dlamch_();

  /*

*

*           Update partial column norms

*

          DO 30 J = I + 1, N

             IF( WORK( J ).NE.ZERO ) THEN

*

*                 NOTE: The following 4 lines follow from the analysis in

*                 Lapack Working Note 176.

*

                TEMP = ABS( A( I, J ) ) / WORK( J )

                TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )

                TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2

                IF( TEMP2 .LE. TOL3Z ) THEN

                   IF( M-I.GT.0 ) THEN

                      WORK( J ) = DNRM2( M-I, A( I+1, J ), 1 )

                      WORK( N+J ) = WORK( J )

                   ELSE

                      WORK( J ) = ZERO

                      WORK( N+J ) = ZERO

                   END IF

                ELSE

                   WORK( J ) = WORK( J )*SQRT( TEMP )

                END IF

             END IF

 30       CONTINUE

  */

  // Some initializations.

  ptr_d  = buff_d;

  ptr_e  = buff_e;

  ptr_wt = buff_wt;

  ptr_A  = buff_A;

  // Main loop.

        } else {

        }

      } else {

      }

    }

  }

}

// ============================================================================

               int m_G, double * buff_G, int ldim_G,

               int pivot_B, int m_B, double * buff_B, int ldim_B,

               int pivot_C, int m_C, double * buff_C, int ldim_C,

               int * buff_p,

               double * buff_d, double * buff_e ) {

//

// It pivots matrix G, pivot vector p, and norms vectors d and e.

// Matrices B and C are optionally pivoted.

//

  double  * ptr_g1, * ptr_g2, * ptr_b1, * ptr_b2, * ptr_c1, * ptr_c2;

  // Swap columns of G, pivots, and norms.

    // Swap full column 0 and column "j_max_col" of G.

    ptr_g1 = & buff_G[ 0 + 0         * ldim_G ];

    // Swap full column 0 and column "j_max_col" of B.

      ptr_b1 = & buff_B[ 0 + 0         * ldim_B ];

    }

    // Swap full column 0 and column "j_max_col" of C.

      ptr_c1 = & buff_C[ 0 + 0         * ldim_C ];

    }

    // Swap element 0 and element "j_max_col" of pivot vector "p".

    // Copy norms of column 0 to column "j_max_col".

  }

}