Dlanst.java

package org.mklab.sdpj.gpack.lapackwrap;

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.gpack.blaswrap.BLAS;


/**
 * @author koga
 * @version $Revision$, 2009/04/24
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
 */
public class Dlanst<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {

  /*  -- LAPACK auxiliary routine (version 3.0) --   
  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
  Courant Institute, Argonne National Lab, and Rice University   
  February 29, 1992   
  
  
  Purpose   
  =======   
  
  DLANST  returns the value of the one norm,  or the Frobenius norm, or   
  the  infinity norm,  or the  element of  largest absolute value  of a   
  real symmetric tridiagonal matrix A.   
  
  Description   
  ===========   
  
  DLANST returns the value   
  
  DLANST = ( max(abs(A(i,j))), NORM = 'M' or 'm'   
           (   
           ( norm1(A),         NORM = '1', 'O' or 'o'   
           (   
           ( normI(A),         NORM = 'I' or 'i'   
           (   
           ( normF(A),         NORM = 'F', 'f', 'E' or 'e'   
  
  where  norm1  denotes the  one norm of a matrix (maximum column sum),   
  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and   
  normF  denotes the  Frobenius norm of a matrix (square root of sum of   
  squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.   
  
  Arguments   
  =========   
  
  NORM    (input) CHARACTER*1   
       Specifies the value to be returned in DLANST as described   
       above.   
  
  N       (input) INTEGER   
       The order of the matrix A.  N >= 0.  When N = 0, DLANST is   
       set to zero.   
  
  D       (input) DOUBLE PRECISION array, dimension (N)   
       The diagonal elements of A.   
  
  RS       (input) DOUBLE PRECISION array, dimension (N-1)   
       The (n-1) sub-diagonal or super-diagonal elements of A.   
  
  =====================================================================   
  */
  /**
   * @param norm norm
   * @param n n
   * @param d d
   * @param e e
   * @return result
   */
  RS execute(String norm, int n, RS[] d, RS[] e) {
    final RS unit = d[0].createUnit();

    int pointer_e = -1;
    int pointer_d = -1;

    RS anorm = null;
    if (n <= 0) {
      anorm = unit.createZero();
    } else if (BLAS.lsame(norm, "M")) { //$NON-NLS-1$

      /* Find max(abs(A(i,j))). */
      RS d1 = d[pointer_d + n];
      anorm = d1.abs();
      for (int i = 1; i <= n - 1; ++i) {
        d1 = d[pointer_d + i];
        anorm = anorm.max(d1.abs());

        d1 = e[pointer_e + i];
        anorm = anorm.max(d1.abs());
        /* L10: */
      }
    } else if (BLAS.lsame(norm, "O") || norm.charAt(0) == '1' || BLAS.lsame(norm, "I")) { //$NON-NLS-1$ //$NON-NLS-2$
      /* Find norm1(A). */
      if (n == 1) {
        anorm = d[pointer_d + 1].abs();
      } else {
        /* Computing MAX */
        RS d3 = (d[pointer_d + 1].abs()).add(e[pointer_e + 1].abs());
        RS d1 = e[pointer_e + n - 1];
        RS d4 = (d1.abs()).add(d[pointer_d + n].abs());
        anorm = d3.max(d4);
        for (int i = 2; i <= n - 1; ++i) {
          /* Computing MAX */
          RS d5 = (d[pointer_d + i].abs()).add(e[pointer_e + i].abs()).add(e[pointer_e + i - 1].abs());
          anorm = anorm.max(d5);
          /* L20: */
        }
      }
    } else if (BLAS.lsame(norm, "F") || BLAS.lsame(norm, "E")) { //$NON-NLS-1$ //$NON-NLS-2$
      /* Find normF(A). */
      RS scale = unit.create(0);
      RS sum = unit.create(1);
      if (n > 1) {
        RS[] eTemp = unit.createArray(e.length - (1 + pointer_e));
        System.arraycopy(e, 1 + pointer_e, eTemp, 0, eTemp.length);
        RS[] ans = Clapack.dlassq(n - 1, eTemp, 1, scale, sum);
        scale = ans[0];
        sum = ans[1];
        sum = sum.multiply(2);
      }
      RS[] dTemp = unit.createArray(d.length - (1 + pointer_d));
      System.arraycopy(d, 1 + pointer_d, dTemp, 0, dTemp.length);
      RS[] ans = Clapack.dlassq(n, dTemp, 1, scale, sum);
      scale = ans[0];
      sum = ans[1];
      anorm = scale.multiply(sum.sqrt());
    }

    RS ret_val = anorm;
    return ret_val;
  }
}