rot.java

JAVA 数学程序库 提供常规的数值计算程序包

      double s, absx, absxy;

      if (xr == 0 && xi==0){
         P.c = 0.;
         P.sr = 1.;
         P.si = 0.;
         P.zr = yr;
         P.zi = yi;
         return;
      }
      s = Math.abs(xr) + Math.abs(xi);
      absx = s*Math.sqrt((xr/s)*(xr/s) + (xi/s)*(xi/s));
      s = Math.abs(s) + Math.abs(yr) + Math.abs(yi);
      absxy = s*Math.sqrt((absx/s)*(absx/s) + (yr/s)*(yr/s) + (yi/s)*(yi/s));
      P.c = absx/absxy;
      xr = xr/absx;
      xi = xi/absx;
      P.sr = -(xr*yr + xi*yi)/absxy;
      P.si = (xi*yr - xr*yi)/absxy;
      P.zr = xr*absxy;
      P.zi = xi*absxy;
   }

/**
    Given a Zmat A, genr returns a plane rotation that on
    postmultiplication into column jj1 and jj2
    annihilates A(ii,jj2).  The element A(ii,jj2) is
    overwirten by zero and the element A(ii,jj1) is
    overwritten by its transformed value.
    @param A    The Zmat (altered)
    @param ii   The index of the row containing the elements
    @param jj1  The column index of the first element
    @param jj2  The column index of the second element (the
                one that is annihilated)
    @return     The rotation
*/

   public static Rot genr(Zmat A, int ii, int jj1, int jj2){

      A.dirty = true;

      int i = ii - A.basex;
      int j1 = jj1 - A.basex;
      int j2 = jj2 - A.basex;

      Rot P = Rot.genr(A.re[i][j1], A.im[i][j1], A.re[i][j2], A.im[i][j2]);
      A.re[i][j1] = P.zr;
      A.im[i][j1] = P.zi;
      A.re[i][j2] = 0;
      A.im[i][j2] = 0;
      return P;
   }

/**
    Given a Zmat A, genr generates a plane rotation that on
    postmultiplication into column jj1 and jj2
    annihilates A(ii,jj2).  The element A(ii,jj2) is
    overwirten by zero and the element A(ii,jj1) is
    overwritten by its transformed value.
    @param A    The Zmat (altered)
    @param ii   The index of the row containing the elements
    @param jj1  The column index of the first element
    @param jj2  The column index of the second element (the
                one that is annihilated)
    @param P    The rotation
*/

   public static void genr(Zmat A, int ii, int jj1, int jj2, Rot P){

      A.dirty = true;

      int i = ii - A.basex;
      int j1 = jj1 - A.basex;
      int j2 = jj2 - A.basex;

      Rot.genr(A.re[i][j1], A.im[i][j1], A.re[i][j2], A.im[i][j2], P);
      A.re[i][j1] = P.zr;
      A.im[i][j1] = P.zi;
      A.re[i][j2] = 0;
      A.im[i][j2] = 0;
   }

/**
   Given  a real 2-vector, genr returns
   a plane rotation such that
<pre>
*      |x y|P = |x y|| c  s||x| = |z 0|
*                    |-s  c||y|
</pre>
   @param x   The first component of the 2-vector
   @param y   The second component of the 2-vector
   @return    The rotation
*/

public static Rot genr(double x, double y){

      Rot P = new Rot();

      P.si = 0.;
      P.zi = 0.;

      double s = Math.abs(x) + Math.abs(y);

      if (s == 0.){
         P.c = 1.;
         P.sr = 0.;
         P.zr = 0.;
         return P;
      }

      P.zr = s*Math.sqrt((x/s)*(x/s) + (y/s)*(y/s));
      P.c = x/P.zr;
      P.sr = -y/P.zr;
      return P;
   }


/**
   Given  a real 2-vector, genr generates
   a plane rotation such that
<pre>
*      |x y|P = |x y|| c  s||x| = |z 0|
*                    |-s  c||y|
</pre>
   @param x   The first component of the 2-vector
   @param y   The second component of the 2-vector
   @param P   The rotation
*/
public static void genr(double x, double y, Rot P){

      P.si = 0.;
      P.zi = 0.;

      double s = Math.abs(x) + Math.abs(y);

      if (s == 0.){
         P.c = 1.;
         P.sr = 0.;
         P.zr = 0.;
         return;
      }

      P.zr = s*Math.sqrt((x/s)*(x/s) + (y/s)*(y/s));
      P.c = x/P.zr;
      P.sr = -y/P.zr;
   }


/**
   Multiplies rows (ii1,jj1:jj2) and (ii2,jj1:jj2)
   of a Zmat (altered) by a plane rotation.
   @param P    The plane rotation
   @param A    The Zmat (altered)
   @param ii1  The row index of the first row.
   @param ii2  The row index of the second row.
   @param jj1  The first index of the range of the rows
   @param jj2  The second index of the range of the rows
*/
   

   public static void pa(Rot P, Zmat A, int ii1, int ii2, int jj1, int jj2 ){

      double t1r, t1i, t2r, t2i;

      A.dirty = true;

      int i1 = ii1 - A.basex;
      int i2 = ii2 - A.basex;
      int j1 = jj1 - A.basex;
      int j2 = jj2 - A.basex;


      for (int j=j1; j<=j2; j++){
         t1r = P.c*A.re[i1][j] + P.sr*A.re[i2][j] - P.si*A.im[i2][j];
         t1i = P.c*A.im[i1][j] + P.sr*A.im[i2][j] + P.si*A.re[i2][j];
         t2r = P.c*A.re[i2][j] - P.sr*A.re[i1][j] - P.si*A.im[i1][j];
         t2i = P.c*A.im[i2][j] - P.sr*A.im[i1][j] + P.si*A.re[i1][j];
         A.re[i1][j] = t1r;
         A.im[i1][j] = t1i;
         A.re[i2][j] = t2r;
         A.im[i2][j] = t2i;
      }  
   }

/**
   Multiplies rows (ii1,jj1:jj2) and (ii2,jj1:jj2)
   of a Zmat (altered) by the conjugate transpose of a plane rotation.
   @param P    The plane rotation
   @param A    The Zmat (altered)
   @param ii1  The row index of the first row.
   @param ii2  The row index of the second row.
   @param jj1  The first index of the range of the rows
   @param jj2  The second index of the range of the rows
*/
   public static void pha(Rot P, Zmat A, int ii1, int ii2, int jj1, int jj2 ){

      double t1r, t1i, t2r, t2i;

      A.dirty = true;

      int i1 = ii1 - A.basex;
      int i2 = ii2 - A.basex;
      int j1 = jj1 - A.basex;
      int j2 = jj2 - A.basex;


      for (int j=j1; j<=j2; j++){
         t1r = P.c*A.re[i1][j] - P.sr*A.re[i2][j] + P.si*A.im[i2][j];
         t1i = P.c*A.im[i1][j] - P.sr*A.im[i2][j] - P.si*A.re[i2][j];
         t2r = P.c*A.re[i2][j] + P.sr*A.re[i1][j] + P.si*A.im[i1][j];
         t2i = P.c*A.im[i2][j] + P.sr*A.im[i1][j] - P.si*A.re[i1][j];
         A.re[i1][j] = t1r;
         A.im[i1][j] = t1i;
         A.re[i2][j] = t2r;
         A.im[i2][j] = t2i;
      }  
   }


/**
   Multiplies columns (ii1:ii2,jj1) and A(ii2:ii2,jj1)
   of a Zmat (altered) by a plane rotation.
   @param A    The Zmat (altered)
   @param P    The rotation
   @param ii1  The first index of the column range
   @param ii2  The second index of the column range
   @param jj1  The index of the first column
   @param jj2  The index of the second column
*/

   public static void ap(Zmat A, Rot P, int ii1, int ii2, int jj1, int jj2 ){

      double t1r, t1i, t2r, t2i;

      A.dirty = true;

      int i1 = ii1 - A.basex;
      int i2 = ii2 - A.basex;
      int j1 = jj1 - A.basex;
      int j2 = jj2 - A.basex;


      for (int i=i1; i<=i2; i++){
         t1r = P.c*A.re[i][j1] - P.sr*A.re[i][j2] - P.si*A.im[i][j2];
         t1i = P.c*A.im[i][j1] - P.sr*A.im[i][j2] + P.si*A.re[i][j2];
         t2r = P.c*A.re[i][j2] + P.sr*A.re[i][j1] - P.si*A.im[i][j1];
         t2i = P.c*A.im[i][j2] + P.sr*A.im[i][j1] + P.si*A.re[i][j1];
         A.re[i][j1] = t1r;
         A.im[i][j1] = t1i;
         A.re[i][j2] = t2r;
         A.im[i][j2] = t2i;
      }  
   }

/**
   Multiplies columns (ii1:ii2,jj1) and A(ii2:ii2,jj1)
   of a Zmat (altered) by the conjugate transpose of plane rotation.
   @param A    The Zmat (altered)
   @param P    The rotation
   @param ii1  The first index of the column range
   @param ii2  The second index of the column range
   @param jj1  The index of the first column
   @param jj2  The index of the second column
*/

   public static void aph(Zmat A, Rot P, int ii1, int ii2, int jj1, int jj2 ){

      double t1r, t1i, t2r, t2i;

      A.dirty = true;

      int i1 = ii1 - A.basex;
      int i2 = ii2 - A.basex;
      int j1 = jj1 - A.basex;
      int j2 = jj2 - A.basex;


      for (int i=i1; i<=i2; i++){
         t1r = P.c*A.re[i][j1] + P.sr*A.re[i][j2] + P.si*A.im[i][j2];
         t1i = P.c*A.im[i][j1] + P.sr*A.im[i][j2] - P.si*A.re[i][j2];
         t2r = P.c*A.re[i][j2] - P.sr*A.re[i][j1] + P.si*A.im[i][j1];
         t2i = P.c*A.im[i][j2] - P.sr*A.im[i][j1] - P.si*A.re[i][j1];
         A.re[i][j1] = t1r;
         A.im[i][j1] = t1i;
         A.re[i][j2] = t2r;
         A.im[i][j2] = t2i;
      }
   }


}