<p>
Inner block mult add operation that takes in account the zeros and on in dataA,
which is the top part of the Y block vector that has the householder vectors.<br>
<br>
C = C + A * B
</p>
9424
public static void multBlockAdd_zerosone( double[] dataA, double []dataB, double []dataC,
                                              int indexA, int indexB, int indexC,
                                              final int heightA, final int widthA, final int widthC) {


        for( int i = 0; i < heightA; i++ ) {
            for( int j = 0; j < widthC; j++ ) {
                double val = i < widthA ? dataB[i*widthC+j+indexB] : 0;

                int end = Math.min(i,widthA);

                for( int k = 0; k < end; k++ ) {
                    val += dataA[i*widthA + k + indexA] * dataB[k*widthC + j + indexB];
                }

                dataC[ i*widthC + j + indexC ] += val;
            }
        }
    }
1262849
public void householderSymmetric( int row , double gamma )
    {
        int startU = row*N;

        // compute v = -gamma*A*u
        for( int i = row+1; i < N; i++ ) {
            double totalReal = 0;
            double totalImag = 0;

            // the lower triangle is not written to so it needs to traverse upwards
            // to get the information.  Reduces the number of matrix writes need
            // improving large matrix performance
            for( int j = row+1; j < i; j++ ) {
                double realA = QT.data[(j*N+i)*2];
                double imagA = -QT.data[(j*N+i)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            for( int j = i; j < N; j++ ) {
                double realA = QT.data[(i*N+j)*2];
                double imagA = QT.data[(i*N+j)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            w[i*2]   = -gamma*totalReal;
            w[i*2+1] = -gamma*totalImag;
        }

        // alpha = -0.5*gamma*u^T*v
        double realAplha = 0;
        double imageAlpha = 0;

        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = -QT.data[(startU+i)*2+1];

            double realV = w[i*2];
            double imagV = w[i*2+1];

            realAplha += realU*realV - imagU*imagV;
            imageAlpha += realU*imagV + imagU*realV;
        }
        realAplha *= -0.5*gamma;
        imageAlpha *= -0.5*gamma;

        // w = v + alpha*u
        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = QT.data[(startU+i)*2+1];

            w[i*2]   += realAplha*realU - imageAlpha*imagU;
            w[i*2+1] += realAplha*imagU + imageAlpha*realU;
        }

        // A = A + w*u^T + u*w^T
        for( int i = row+1; i < N; i++ ) {

            double realWW = w[i*2];
            double imagWW = w[i*2+1];

            double realUU = QT.data[(startU+i)*2];
            double imagUU = QT.data[(startU+i)*2+1];

            int indA = (i*N+i)*2;
            for( int j = i; j < N; j++ ) {
                // only write to the upper portion of the matrix
                // this reduces the number of cache misses

                double realU = QT.data[(startU+j)*2];
                double imagU = -QT.data[(startU+j)*2+1];

                double realW = w[j*2];
                double imagW = -w[j*2+1];

                QT.data[indA++] += realWW*realU - imagWW*imagU + realW*realUU - imagW*imagUU;
                QT.data[indA++] += realWW*imagU + imagWW*realU + realW*imagUU + imagW*realUU;
            }
        }
    }
1256470
public static void checkVectorAddition() {
        DoubleMatrix x = new DoubleMatrix(3, 1, 1.0, 2.0, 3.0);
        DoubleMatrix y = new DoubleMatrix(3, 1, 4.0, 5.0, 6.0);
        DoubleMatrix z = new DoubleMatrix(3, 1, 5.0, 7.0, 9.0);

        check("checking vector addition", x.add(y).equals(z));
    }
1248419
public Matrix plus(Matrix B) {
        Matrix A = this;
        if (B.M != A.M || B.N != A.N) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }
        Matrix C = new Matrix(M, N);
        for (int i = 0; i < M; i++) {
            for (int j = 0; j < N; j++) {
                C.data[i][j] = A.data[i][j] + B.data[i][j];
            }
        }
        return C;
    }
73307
public static void addCol(Matrix A, int j, double c)
    {
        addCol(A, j, 0, A.rows(), c);
    }
1303325
public static void symmRankNMinus_U(int blockLength ,
                                        DSubmatrixD1 A , DSubmatrixD1 B )
    {

        int heightB = B.row1-B.row0;
        if( heightB > blockLength )
            throw new IllegalArgumentException("Height of B cannot be greater than the block length");

        int N = B.col1-B.col0;

        if( A.col1-A.col0 != N )
            throw new IllegalArgumentException("A does not have the expected number of columns based on B's width");
        if( A.row1-A.row0 != N )
            throw new IllegalArgumentException("A does not have the expected number of rows based on B's width");


        for( int i = B.col0; i < B.col1; i += blockLength ) {

            int indexB_i = B.row0*B.original.numCols + i*heightB;
            int widthB_i = Math.min(blockLength,B.col1-i);

            int rowA = i-B.col0+A.row0;
            int heightA = Math.min( blockLength , A.row1 - rowA);

            for( int j = i; j < B.col1; j += blockLength ) {

                int widthB_j = Math.min(blockLength,B.col1-j);

                int indexA = rowA * A.original.numCols + (j-B.col0+A.col0)*heightA;
                int indexB_j = B.row0*B.original.numCols + j*heightB;

                if( i == j ) {
                    // only the upper portion of this block needs to be modified since it is along a diagonal
                    multTransABlockMinus_U( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,heightB,widthB_i,widthB_j);
                } else {
                    multTransABlockMinus( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,heightB,widthB_i,widthB_j);
                }
            }
        }
    }
747566
public Matrix add(Matrix B)
    {
        Matrix toReturn = getThisSideMatrix(B);
        toReturn.mutableAdd(1.0, B);
        return toReturn;
    }
1258424
public Vector4d add(double x, double y, double z, double w) {
        return add(x, y, z, w, thisOrNew());
    }
1344386
public Vector add(Vector b) {
        if ((b == null) || (size != b.size)) return null;

        int i;
        Vector result = new Vector(size);
        for (i = 0; i < size; i++)
            result.vector[i] = vector[i] + b.vector[i];
        return result;
    }
154021
@Override
    public QR qr() {
        // Initialize.
        int m = nrows();
        int n = ncols();
        double[] rDiagonal = new double[n];

        // Main loop.
        for (int k = 0; k < n; k++) {
            // Compute 2-norm of k-th column without under/overflow.
            double nrm = 0.0;
            for (int i = k; i < m; i++) {
                nrm = Math.hypot(nrm, get(i, k));
            }

            if (nrm != 0.0) {
                // Form k-th Householder vector.
                if (get(k, k) < 0) {
                    nrm = -nrm;
                }
                for (int i = k; i < m; i++) {
                    div(i, k, nrm);
                }
                add(k, k, 1.0);

                // Apply transformation to remaining columns.
                for (int j = k + 1; j < n; j++) {
                    double s = 0.0;
                    for (int i = k; i < m; i++) {
                        s += get(i, k) * get(i, j);
                    }
                    s = -s / get(k, k);
                    for (int i = k; i < m; i++) {
                        add(i, j, s * get(i, k));
                    }
                }
            }
            rDiagonal[k] = -nrm;
        }

        boolean singular = false;
        for (int j = 0; j < rDiagonal.length; j++) {
            if (rDiagonal[j] == 0) {
                singular = true;
                break;
            }
        }

        return new QR(this, rDiagonal, singular);
    }
994844
private void addToJacobian(DMatrix tripplet, int col , double a[], double b[], DMatrixRMaj R ) {
		set(tripplet,jacRowX,col+0,a[0]*R.data[0] + a[1]*R.data[3] + a[2]*R.data[6]);
		set(tripplet,jacRowX,col+1,a[0]*R.data[1] + a[1]*R.data[4] + a[2]*R.data[7]);
		set(tripplet,jacRowX,col+2,a[0]*R.data[2] + a[1]*R.data[5] + a[2]*R.data[8]);

		set(tripplet,jacRowY,col+0,b[0]*R.data[0] + b[1]*R.data[3] + b[2]*R.data[6]);
		set(tripplet,jacRowY,col+1,b[0]*R.data[1] + b[1]*R.data[4] + b[2]*R.data[7]);
		set(tripplet,jacRowY,col+2,b[0]*R.data[2] + b[1]*R.data[5] + b[2]*R.data[8]);
	}
44515
private final void MultiplyAv(int n, double[] v, double[] Av){
        for (int i=0; i<n; i++){
            Av[i] = 0;
            for (int j=0; j<n; j++) Av[i] += A(i,j)*v[j];
        }
    }

     transposed */
    private final void MultiplyAtv(int n, double[] v, double[] Atv){
        for (int i=0;i<n;i++){
            Atv[i] = 0;
            for (int j=0; j<n; j++) Atv[i] += A(j,i)*v[j];
        }
    }

     and then by matrix A transposed */
    private final void MultiplyAtAv(int n, double[] v, double[] AtAv){
        double[] u = new double[n];
        MultiplyAv(n,v,u);
        MultiplyAtv(n,u,AtAv);
    }
}
417270
public void mul(double b, IVector result) {
        if (result.size != size) result.reshape(size);

        int i;
        for (i = 0; i < size; i++) {
            result.realvector[i] = realvector[i] * b;
            result.imagvector[i] = imagvector[i] * b;
        }
    }
1294034
public void mul(IMatrix b, IMatrix result) {
        if ((b == null) || (columns != b.rows)) return;

        if ((result.rows != rows) || (result.columns != b.columns)) result.reshape(rows, b.columns);
        int i, j, k;
        double realsum, imagsum;
        for (i = 0; i < rows; i++)
            for (k = 0; k < b.columns; k++) {
                realsum = 0;
                imagsum = 0;
                for (j = 0; j < columns; j++) {
                    realsum += realmatrix[i][j] * b.realmatrix[j][k] - imagmatrix[i][j] * b.imagmatrix[j][k];
                    imagsum += realmatrix[i][j] * b.imagmatrix[j][k] + imagmatrix[i][j] * b.realmatrix[j][k];
                }
                result.realmatrix[i][k] = realsum;
                result.imagmatrix[i][k] = imagsum;
            }
    }
1240852
public static void byRow(Matrix m) {
        for (int i = 0; i < m.rows(); ++i) {
            double rowSum = 0;
            for (int j = 0; j < m.columns(); ++j)
                rowSum += m.get(i,j);

            // Skip row whose sume is zero.
            if (rowSum == 0) 
                continue;
            for (int j = 0; j < m.columns(); ++j)
                m.set(i,j, m.get(i,j) / rowSum);
        }
    }
1337353
void incr1( int b, double y, double yy ) {
    Utils.AtomicDoubleArray.add(_sums,b,y);
    Utils.AtomicDoubleArray.add(_ssqs,b,yy);
  }
439096
public static boolean bidiagOuterBlocks( final int blockLength ,
                                             final DSubmatrixD1 A ,
                                             final double gammasU[],
                                             final double gammasV[])
    {
//        System.out.println("---------- Orig");
//        A.original.print();

        int width = Math.min(blockLength,A.col1-A.col0);
        int height = Math.min(blockLength,A.row1-A.row0);

        int min = Math.min(width,height);

        for( int i = 0; i < min; i++ ) {
            //--- Apply reflector to the column

            // compute the householder vector
            if (!computeHouseHolderCol(blockLength, A, gammasU, i))
                return false;

            // apply to rest of the columns in the column block
            rank1UpdateMultR_Col(blockLength,A,i,gammasU[A.col0+i]);

            // apply to the top row block
            rank1UpdateMultR_TopRow(blockLength,A,i,gammasU[A.col0+i]);

            System.out.println("After column stuff");
            A.original.print();

            //-- Apply reflector to the row
            if(!computeHouseHolderRow(blockLength,A,gammasV,i))
                return false;
            
            // apply to rest of the rows in the row block
            rank1UpdateMultL_Row(blockLength,A,i,i+1,gammasV[A.row0+i]);

            System.out.println("After update row");
            A.original.print();

            // apply to the left column block
            // TODO THIS WON'T WORK!!!!!!!!!!!!!
            // Needs the whole matrix to have been updated by the left reflector to compute the correct solution
//            rank1UpdateMultL_LeftCol(blockLength,A,i,i+1,gammasV[A.row0+i]);

            System.out.println("After row stuff");
            A.original.print();
        }

        return true;
    }