/*## 
 ##  class NMatrix
 ##   
 ##  This class contains all the BLAS (basic linear algebra subroutines)
 ##   and other low-level functions common to all matrices.  
 ##  NMatrix contains no factoring or solving capabilities.
 ##
 ##*/

package edu.rice.linpack.Matrix.NMatrix;
import java.lang.*;
import edu.rice.linpack.LNumber.*;
import edu.rice.linpack.util.*;

public abstract class NMatrix {

  protected int rows, cols;
  protected int[] pivot;
  protected LNumber[][] Mat;
  
  /*##
   ##  These are the default methods for the factoring routines common
   ##   to most subclasses.  If a routine is not supported by a 
   ##   subclass due to the storage format, an UnsupportedMethodException 
   ##   is thrown.
   ##*/

  public void factor() 
       throws SingularMatrixException, UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public LNumber condition()
       throws SingularMatrixException, UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public LNumber condition(LNumber[] Z)
       throws SingularMatrixException, UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public void solve(LNumber[] B, int i)
       throws UnsupportedMethodException,   
	 UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public void solve(LNumber[] B)
       throws SingularMatrixException, UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public void inverse() 
       throws UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public LNumber[] determ() 
       throws UnsupportedMethodException
  {
    throw new UnsupportedMethodException();
  }
  public void Factor_and_Solve(LNumber[] B) 
       throws SingularMatrixException, UnsupportedMethodException
  {
    this.factor();
    this.solve(B);
  }

  /*##  
   ##   These are used as the constructors for several subclasses
   ##*/
  
  void cons() {
    pivot = new int[cols];
  }
  void cons(int i, int j) {
    Mat = new LNumber[i][j];
    rows = i;
    cols = j;
    pivot = new int[cols];
  }
  void cons(NMatrix M) {
    rows = M.numofRows();
    cols = M.numofCols();
    Mat = new LNumber[rows][cols];
    for(int i=0;i<rows;i++)
      for(int j=0;j<cols;j++) 
	Mat[i][j] = (M.getElem(i,j)).Clone();
    pivot = new int[cols];
  }    
  void cons(LNumber[][] F) {
    rows = F.length;
    cols = (F[0]).length;
    Mat = new LNumber[rows][cols];
    Mat = F;
    pivot = new int[cols];
  }

  public int numofRows() {
    return rows;
  }
  public int numofCols() {
    return cols;
  }
  public int[] getPivots() {
    return pivot;
  }  

  /*##
   ##  These routines return a copy of the desired object beginning at
   ##   the element (i,q)  
   ##*/

  public LNumber getElem(int i, int q) {
    return Mat[i][q].Clone();
  }
  public LNumber[] getRow(int i) {
    return this.getRow(i,0);
  }
  public LNumber[] getRow(int i, int q) {
    LNumber[] F = new LNumber[cols-q];
    for(int j=0;j<(cols-q);j++) 
      F[j] = Mat[i][q+j].Clone();
    return F;
  }
  public LNumber[] getColumn(int i) {
    return this.getColumn(i,0);
  }
  public LNumber[] getColumn(int i, int q) {
    LNumber[] F = new LNumber[rows-q];
    for(int j=0; j<(rows-q);j++) {
      F[j] = Mat[j+q][i].Clone();
    }
    return F;
  }
  public LNumber newOfType() {
    return Mat[0][0].Clone();
  }

  /*##
   ##  this.setColumn changes column c to F  
   ##*/

  public void setColumn(int c,LNumber[] F) {
    for(int r=0;r<rows;r++) 
      Mat[r][c] = F[r].Clone();
  }

  /*## 
   ## this.setElem changes element r,c to value
   ##*/

  public void setElem(int r, int c, LNumber value) {
    Mat[r][c] = value.Clone();
  }

  public void Print() {
    System.out.println(); 
    for(int i=0;i<rows;i++) {
      for(int j=0;j<cols;j++) {
	if(Mat[i][j] != null) {
	  Mat[i][j].Print();
	  System.out.print("   ");
	}
	else {
	  System.out.print("0.0    ");
	}
      }
      System.out.println();
    }
    System.out.println();
  }

  /*##
   ##  this.swapElem swaps two elements in this NMatrix  
   ##*/

  public void swapElem(int rx, int cx, int ry, int cy) {
    LNumber T = Mat[rx][cx].Clone();
    Mat[rx][cx] = Mat[ry][cy];
    Mat[ry][cy] = T;
  }
  

  /*#######################*/
  /*##   BLAS Routines   ##*/

  /*##
   ##  For all BLAS: 
   ##
   ##  On Input:
   ##
   ##     n         int
   ##               The number of elements to be effected. 
   ##
   ##     incx      int
   ##               The amount to increase the column number for 
   ##               this.Mat on each iteration.
   ##
   ##     r,rx      int 
   ##               The initial row in this.Mat.
   ##
   ##     c,cx      int
   ##               The initial column in this.Mat.
   ## 
   ##     B         LNumber[], NMatrix, LNumber[][]
   ##               The other object used in some calculations.
   ##     
   ##     incy      int
   ##               The amount to increase the column number for B.
   ##
   ##     ry        int
   ##               The initial row in B if B is an NMatrix or LNumber[][].
   ## 
   ##     cy        int
   ##               The initial column in B.
   ##
   ##     Da        LNumber
   ##               A LNumber used as a scalar.
   ##
   ##     Tg        Trig
   ##               Holds the sine and cosine values formed in a Givens
   ##               rotation and used in a rotation.
   ##*/

  /*##
   ##  asum - sum of absolute values in a vector
   ##  Returns the sum as a LNumber  
   ##*/

  public LNumber asum(int n, int incx, int r, int c) {
    
    LNumber dasum = Mat[r][c].abs();
    if(n > 0) {
      if(incx != 1) {
	int ix = r+incx;
	for(int i=1;i<n;i++,ix+=incx) 
	  dasum.addTo(Mat[ix][c].abs());
      }
      else {
	int m = (n-1)%6;
	if(m != 0) {
	  for(int i=1;i<=m;i++) 
	    dasum.addTo(Mat[i+r][c].abs());
	}
	if(n >= 6) {
	  for(int i=m+1;i<n;i+=6) {
	    
	    int ir = i+r;
	    dasum.addTo((((((Mat[ir][c].abs()).add(Mat[ir+1][c].abs()))
			   .add(Mat[ir+2][c].abs())).add(Mat[ir+3][c].abs()))
			 .add(Mat[ir+4][c].abs())).add(Mat[ir+5][c].abs()));
	  }
	}
      }
    }
    else 
      dasum.setZero();
    return dasum;
  }

  /*##
   ##  axpy - a const times a vector plus a vector 
   ##  B is changed in the loop and holds the desired value on return 
   ##*/

  public void axpy(int n, LNumber Da, int incx, int rx, int cx, 
		   NMatrix B, int incy, int ry, int cy) {
    
    if(n > 0 && !Da.equals(0)) {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	for(int i=0;i<n;i++,ix+=incx,iy+=incy) 
	  B.setElem(iy+ry,cy,
		    (B.getElem(iy+ry,cy)).add(Da.mult(Mat[ix+rx][cx])));
      }
      else {
	int m = n%4;
        if(m != 0) {
          for(int i=0;i<m;i++)
            B.setElem(i+ry,cy,
		      ((B.getElem(i+ry,cy)).add(Da.mult(Mat[i+rx][cx]))));
	}
        if(n >= 4) {
          for(int i=m;i<n;i+=4) {
            int iry = i+ry;
            int irx = i+rx;
            B.setElem(iry,cy,
		      (B.getElem(iry,cy)).add(Da.mult(Mat[irx][cx])));
	    B.setElem(iry+1,cy,
		      (B.getElem(iry+1,cy)).add(Da.mult(Mat[irx+1][cx])));
	    B.setElem(iry+2,cy,
		      (B.getElem(iry+2,cy)).add(Da.mult(Mat[irx+2][cx])));
	    B.setElem(iry+3,cy,
		      (B.getElem(iry+3,cy)).add(Da.mult(Mat[irx+3][cx])));
          }
        }
      }
    }
  }
  public void axpy(int n, LNumber Da, int incx, int rx, int cx,
		   LNumber[] B, int incy, int ry) {
    
    if(n > 0 && !Da.equals(0)) {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	for(int i=0;i<n;i++,ix+=incx,iy+=incy) 
	  B[ry+iy].addTo(Da.mult(Mat[ix+rx][cx]));
      }
      else {
	int m = n%4;
        if(m != 0) {
          for(int i=0;i<m;i++)
            B[i+ry].addTo(Da.mult(Mat[i+rx][cx]));
        }
        if(n >= 4) {
          for(int i=m;i<n;i+=4) {
            int iry = i+ry;
            int irx = i+rx;
            B[iry].addTo(Da.mult(Mat[irx][cx]));
            B[iry+1].addTo(Da.mult(Mat[irx+1][cx]));
            B[iry+2].addTo(Da.mult(Mat[irx+2][cx]));
            B[iry+3].addTo(Da.mult(Mat[irx+3][cx]));
          }
        }
      }
    }
  }
  
  /*##
   ##  This version takes an array of LNumbers and modifies the calling object.
   ##  this.Mat is changed in the loop and holds the desired value on return.
   ##*/

  public void axpy(int n, LNumber Da, LNumber[] B, int incy, int ry,
		   int incx, int rx, int cx) {
    if(n > 0 && !Da.equals(0)) {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	for(int i=0;i<n;i++,ix+=incx,iy+=incy) 
	  Mat[rx+ix][cx].addTo(Da.mult(B[ry+iy]));
      }
      else {
	int m = n%4;
        if(m != 0) {
          for(int i=0;i<m;i++)
            Mat[i+rx][cx].addTo(Da.mult(B[i+ry]));
        }
        if(n >= 4) {
          for(int i=m;i<n;i+=4) {
            int iry = i+ry;
            int irx = i+rx;
            Mat[irx][cx].addTo(Da.mult(B[iry]));
            Mat[irx+1][cx].addTo(Da.mult(B[iry+1]));
            Mat[irx+2][cx].addTo(Da.mult(B[iry+2]));
            Mat[irx+3][cx].addTo(Da.mult(B[iry+3]));
          }
        }
      }
    }
  }


  /*##
   ##  copy - copies one vector to another 
   ##  B is changed and holds the desired array 
   ##*/

  public void copy(int n, int incx, int rx, int cx, 
		   LNumber[] B, int incy) {
    if(n > 0) {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	
	for(int i=0;i<n;i++,ix+=incx, iy+= incy) 
	  B[iy] = Mat[ix+rx][cx].Clone();
      }
      else {
	int m = n%7;
	if(m != 0) {
          for(int i=0;i<m;i++) 
            B[i] = Mat[i+rx][cx].Clone();
        }
        if(n >= 7) {
          for(int i=m;i<n;i+=7) {
            int irx = i+rx;
            B[i] = Mat[irx][cx].Clone();
            B[i+1] = Mat[irx+1][cx].Clone();
            B[i+2] = Mat[irx+2][cx].Clone();
            B[i+3] = Mat[irx+3][cx].Clone();
            B[i+4] = Mat[irx+4][cx].Clone();
            B[i+5] = Mat[irx+5][cx].Clone();
            B[i+6] = Mat[irx+6][cx].Clone();
          }
        }
      }
    }
  }
  public void copy(int n, int incx, int rx, int cx, 
		   NMatrix B, int incy, int ry, int cy) {
    
    if(n > 0) {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	
	for(int i=0;i<n;i++,ix+=incx, iy+= incy) 
	  B.setElem(iy+ry,cy,Mat[ix+rx][cx].Clone());
      }
      else {
	        int m = n%7;
        if(m != 0) {
          for(int i=0;i<m;i++) 
	    B.setElem(i+ry,cy,Mat[i+rx][cx].Clone());
	}
        if(n >= 7) {
          for(int i=m;i<n;i+=7) {
            int irx = i+rx;
            int iry = i+ry;
	    B.setElem(iry,cy,Mat[irx][cx].Clone());
	    B.setElem(iry+1,cy,Mat[irx+1][cx].Clone());
	    B.setElem(iry+2,cy,Mat[irx+2][cx].Clone());
	    B.setElem(iry+3,cy,Mat[irx+3][cx].Clone());
	    B.setElem(iry+4,cy,Mat[irx+4][cx].Clone());
	    B.setElem(iry+5,cy,Mat[irx+5][cx].Clone());
	    B.setElem(iry+6,cy,Mat[irx+6][cx].Clone());
          }
	}
      }
    }
  }
  

  /*##
   ##  dot - Dot product between two vectors 
   ##  An LNumber is returned  
   ##*/

  public LNumber dot(int n, int incx, int rx, int cx,
		   NMatrix B, int incy, int ry, int cy) {
    
    LNumber N = Mat[rx][cx].Clone();
    N.setZero();
    
    if(n > 0){
      if(incy != 1 || incx != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	
	// CHANGING THIS - ADDING CONJ()

	for(int i=0; i<n;i++,ix+= incx,iy += incy) 
	  N.addTo((Mat[rx+ix][cx].conj()).mult(B.getElem(ry+iy,cy)));
      }
      else {
	int m = n%5;
        if(m != 0) {
          for(int i=0;i<m;i++) {

	// CHANGING THIS - ADDING CONJ()

	    N.addTo((Mat[rx+i][cx].conj()).mult(B.getElem(ry+i,cy)));
          }
        }
        if(n >= 5) {
          for(int i=m;i<n;i+=5) {
            int irx = i+rx;
            int iry = i+ry;

	// CHANGING THIS - ADDING CONJ()

            N.addTo((((((Mat[irx][cx].conj()).mult(B.getElem(iry,cy)))
		       .add((Mat[irx+1][cx].conj()).mult(B.getElem(iry+1,cy))))
		      .add((Mat[irx+2][cx].conj()).mult(B.getElem(iry+2,cy))))
		     .add((Mat[irx+3][cx].conj()).mult(B.getElem(iry+3,cy))))
		    .add((Mat[irx+4][cx].conj()).mult(B.getElem(iry+4,cy))));
          }
        }
      }
    }
    return N;
  }
  public LNumber dot(int n, int incx, int rx, int cx,
		   LNumber[] B, int incy, int ry) {

    LNumber N = Mat[rx][cx].Clone();
    N.setZero();
    
    if(n > 0){
      if(incy != 1 || incx != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	
	for(int i=0; i<n;i++,ix+= incx,iy += incy) 

	// CHANGING THIS - ADDING CONJ()

	  N.addTo((Mat[rx+ix][cx].conj()).mult(B[ry+iy]));
      }
      else {
	int m = n%5;
        if(m != 0) {
          for(int i=0;i<m;i++) {

	// CHANGING THIS - ADDING CONJ()

            N.addTo((Mat[i+rx][cx].conj()).mult(B[ry+i]));
          }
        }
        if(n >= 5) {
          for(int i=m;i<n;i+=5) {
            int irx = i+rx;
            int iry = i+ry;

	// CHANGING THIS - ADDING CONJ()

            N.addTo((((((Mat[irx][cx].conj()).mult(B[iry])) 
		       .add((Mat[irx+1][cx].conj()).mult(B[iry+1])))
		      .add((Mat[irx+2][cx].conj()).mult(B[iry+2])))
		     .add((Mat[irx+3][cx].conj()).mult(B[iry+3])))
		    .add((Mat[irx+4][cx].conj()).mult(B[iry+4])));
          }
        }
      }
    }
    return N;
  }


  /*##
   ##  scal - multiply a vector by a scalar 
   ##  The calling object is changed  
   ##*/

  public void scal(int n, LNumber Da, int incx, int r, int c) {

    if(n > 0 && !Da.equals(1)) {
      if(incx != 1) {
	int ix = 0;
	for (int i=0;i<n;i++,ix += incx)
	  Mat[ix+r][c].multTo(Da);
      }
      else {
	int m = n%5;
        if(m != 0) {
          for(int i=0;i<m;i++) {
            Mat[i+r][c].multTo(Da);
          }
        }
        if(n >= 5) {
          for(int i=m;i<n;i+=5) {
            int ir = i+r;
            Mat[ir][c].multTo(Da); 
            Mat[ir+1][c].multTo(Da);
            Mat[ir+2][c].multTo(Da); 
            Mat[ir+3][c].multTo(Da); 
            Mat[ir+4][c].multTo(Da);
          }
        }
      }
    }      
  }      
  
  public void scal(int n, int scale, int incx, int r, int c) {

    if(n > 0 && (scale != 1)) {
      if(incx != 1) {
	int ix = 0;
	for (int i=0;i<n;i++,ix += incx)
	  Mat[ix+r][c].multTo(scale);
      }
      else {
	int m = n%5;
        if(m != 0) {
          for(int i=0;i<m;i++) {
            Mat[i+r][c].multTo(scale);
          }
        }
        if(n >= 5) {
          for(int i=m;i<n;i+=5) {
            int ir = i+r;
            Mat[ir][c].multTo(scale); 
            Mat[ir+1][c].multTo(scale);
            Mat[ir+2][c].multTo(scale); 
            Mat[ir+3][c].multTo(scale); 
            Mat[ir+4][c].multTo(scale);
          }
        }
      }
    }      
  }      
  

  /*##
   ##  swap - interchanges two vectors 
   ##*/

  public void swap(int n, int incx, int rx, int cx, 
		   NMatrix B, int incy, int ry, int cy) {
    
    if(n > 0) {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	LNumber N;
	for(int i=0;i<n;i++) {
	  N = Mat[ix+rx][cx].Clone();
	  Mat[ix+rx][cx] = B.getElem(iy+ry,cy);
	  B.setElem(iy+ry,cy,N);
	  ix += incx;
	  iy += incy;
	}
      }
      else {
	int m = n%3;
        if(m != 0) {
          LNumber temp;
          for(int i=0;i<m;i++) {
            temp = Mat[i+rx][cx].Clone();
            Mat[i+rx][cx] = B.getElem(i+ry,cy);
	    B.setElem(i+ry,cy,temp);
          }
        }
        if(n >= 3) {
          LNumber temp;
          for(int i=m;i<n;i+=3) {
            int irx = i+rx;
            int iry = i+ry;
            temp = Mat[irx][cx].Clone();
            Mat[irx][cx] = B.getElem(iry,cy);
	    B.setElem(iry,cy,temp);
            temp = Mat[irx+1][cx].Clone();
            Mat[irx+1][cx] = B.getElem(iry+1,cy);
	    B.setElem(iry+1,cy,temp);
            temp = Mat[irx+2][cx].Clone();
            Mat[irx+2][cx] = B.getElem(iry+2,cy);
	    B.setElem(iry+2,cy,temp);
          }
        }
      }
    }
  }

  /*##
   ##  i_amax - finds the index of element having the max absolute value 
   ##  Returns the index with zero as first element 
   ##*/

  public int i_amax(int n, int incx, int r, int c) {
  
    int max = 0;
    LNumber Dmax;
    
    if(n < 0) 
      return 0;  
    
    if(n > 0) {
      if(incx != 1) {
	int ix = incx;
	Dmax = Mat[r][c].abs();
	for(int i=1;i<n;i++) {
	  if(!((Mat[ix+r][c].abs()).lessOrEqual(Dmax))) {
	    max = ix;
	    Dmax = Mat[ix+r][c].abs();
	  }
	  ix += incx;
	}
      } 
      else {
	Dmax = Mat[r][c].abs();
        for(int i=1;i<n;i++) {
          if(!((Mat[r+i][c].abs()).lessOrEqual(Dmax))) {
            max = i;
            Dmax = Mat[r+i][c].abs();
          }
        }
      }
    }
    return max;
  }
  
  /*##
   ##  rotg - construct a Givens plane rotation 
   ##  Tger contains the computed sin and cos values on return
   ##*/
  
  public void rotg(int r, int c, Trig Tger) {
    
    NTrig Tg = (NTrig)Tger;

    LNumber Roe = Tg.getB();
    LNumber Da = this.Mat[r][c];
    LNumber R;

    if((Da.abs()).greaterThan((Tg.getB()).abs())) 
      Roe = Da;
    LNumber Scale = (Da.abs()).add((Tg.getB().abs()));
    if(Scale.equals(0)) {
      Tg.setCosOne();
      Tg.setSinZero();
      R = Scale.Clone();
      R.setZero();
    }
    else {
      R = Scale.mult((((Da.div(Scale)).square())
		      .add(((Tg.getB()).div(Scale).square()))).sqrt());
      if(!Roe.greaterOrEqual(0)) 
	R.negateTo();
      Tg.setCos(Da.div(R));
      Tg.setSin(Tg.getB().div(R));
    }
    LNumber Z;
    if((Da.abs()).greaterThan((Tg.getB()).abs())) 
      Z = (Tg.getSin()).Clone();
    else if(!(Tg.getCos()).equals(0)) 
      Z = (Tg.getCos()).inv();
    else {
      Z = Roe.Clone();
      Z.setOne();
    }
    Mat[r][c] = R.Clone();
    Tg.setB(Z);
  }
  
  /*##
   ##  rot applies a plane rotation from the parameters figured from Givens 
   ##*/

  public void rot(int n, int incx, int rx, int cx, NMatrix B, 
		  int incy, int ry, int cy, Trig Tger) {
    
    NTrig Tg = (NTrig) Tger;
    if(n <= 0) 
      return;
    else {
      if(incx != 1 || incy != 1) {
	int ix = 0;
	int iy = 0;
	if(incx < 0) ix = (-n+1)*incx;
	if(incy < 0) iy = (-n+1)*incy;
	
	LNumber temp;      
	for(int i=0;i<n;i++) {
	  temp = (Mat[ix+rx][cx].mult(Tg.getCos()))
	    .add((B.getElem(iy+ry,cy)).mult(Tg.getSin()));
	  B.setElem(iy+ry,cy,((B.getElem(iy+ry,cy)).mult(Tg.getCos()))
		    .sub(Mat[ix+rx][cx].mult(Tg.getSin())));
	  Mat[ix+rx][cx] = temp.Clone();
	  ix += incx;
	  iy += incy;
	}
      }
      else {
	LNumber temp;
        for(int i=0;i<n;i++) {
          temp = (Mat[i+rx][cx].mult(Tg.getCos()))
	    .add((B.getElem(i+ry,cy)).mult(Tg.getSin()));
      	  B.setElem(i+ry,cy,((B.getElem(i+ry,cy)).mult(Tg.getCos()))
		    .sub(Mat[i+rx][cx].mult(Tg.getSin())));
          Mat[i+rx][cx] = temp.Clone();
        }
      }
    }
  }

  /*##
   ##  Finds the 2-norm of a column in an NMatrix
   ##  Returns the value as an LNumber
   ##*/

  public LNumber nrm2(int n, int incx, int r, int c) {
    
    int ix;
    LNumber absxi, scale, ssq;
    scale = Mat[r][c].Clone();
    scale.setZero();

    if(n < 1 || incx < 1) 
      return scale;
    else if( n == 1) 
      return Mat[r][c].abs();
    
    ssq = scale.Clone();
    ssq.setOne();
    
    for(ix = 0;ix < n*incx;ix+=incx) 
      Mat[ix+r][c].ssq(ssq,scale);

    return (scale.mult(ssq.sqrt()));
  }
}