Performs a standard bidiagonal decomposition just on the outer blocks of the provided matrix

@param blockLength
@param A
@param gammasU
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;
    }
18220
private static void calculateThresholdForBlock(byte[] luminances,
                                                 int subWidth,
                                                 int subHeight,
                                                 int width,
                                                 int height,
                                                 int[][] blackPoints,
                                                 BitMatrix matrix) {
    int maxYOffset = height - BLOCK_SIZE;
    int maxXOffset = width - BLOCK_SIZE;
    for (int y = 0; y < subHeight; y++) {
      int yoffset = y << BLOCK_SIZE_POWER;
      if (yoffset > maxYOffset) {
        yoffset = maxYOffset;
      }
      int top = cap(y, 2, subHeight - 3);
      for (int x = 0; x < subWidth; x++) {
        int xoffset = x << BLOCK_SIZE_POWER;
        if (xoffset > maxXOffset) {
          xoffset = maxXOffset;
        }
        int left = cap(x, 2, subWidth - 3);
        int sum = 0;
        for (int z = -2; z <= 2; z++) {
          int[] blackRow = blackPoints[top + z];
          sum += blackRow[left - 2] + blackRow[left - 1] + blackRow[left] + blackRow[left + 1] + blackRow[left + 2];
        }
        int average = sum / 25;
        thresholdBlock(luminances, xoffset, yoffset, average, width, matrix);
      }
    }
  }
1258029
public static void blockMultPlusTransA( final double[] dataA, final double []dataB, final double []dataC,
                                     int indexA, int indexB, int indexC,
                                     final int heightA, final int widthA, final int widthC) {
//        for( int i = 0; i < widthA; i++ ) {
//            for( int k = 0; k < heightA; k++ ) {
//                double valA = dataA[k*widthA + i + indexA];
//                for( int j = 0; j < widthC; j++ ) {
//                    dataC[ i*widthC + j + indexC ] += valA * dataB[k*widthC + j + indexB];
//                }
//            }
//        }

        int rowC = indexC;
        for( int i = 0; i < widthA; i++ , rowC += widthC) {
            int colA = i + indexA;
            int endA = colA + widthA*heightA;
            int b = indexB;

            // for( int k = 0; k < heightA; k++ ) {
            while(colA != endA ) {
                double valA = dataA[colA];

                int c = rowC;
                final int endB = b + widthC;

                //for( int j = 0; j < widthC; j++ ) {
                while( b != endB ) {
                    dataC[ c++ ] += valA * dataB[b++];
                }
                colA += widthA;
            }
        }
    }
1042886
private static void thresholdBlock(byte[] luminances,
                                     int xoffset,
                                     int yoffset,
                                     int threshold,
                                     int stride,
                                     BitMatrix matrix) {
    for (int y = 0, offset = yoffset * stride + xoffset; y < BLOCK_SIZE; y++, offset += stride) {
      for (int x = 0; x < BLOCK_SIZE; x++) {
        // Comparison needs to be <= so that black == 0 pixels are black even if the threshold is 0.
        if ((luminances[offset + x] & 0xFF) <= threshold) {
          matrix.set(xoffset + x, yoffset + y);
        }
      }
    }
  }
975232
void init( DMatrixSparseCSC A ) {
        this.A = A;
        this.m = A.numRows;
        this.n = A.numCols;

        this.next = 0;
        this.head = m;
        this.tail = m + n;
        this.nque = m + 2*n;

        if( parent.length < n || leftmost.length < m) {
            parent = new int[n];
            post = new int[n];
            pinv = new int[m+n];
            countsR = new int[n];
            leftmost = new int[m];
        }
        gwork.reshape(m+3*n);
    }
366364
public void applyQ(DMatrixRBlock B , boolean isIdentity ) {
        int minDimen = Math.min(dataA.numCols,dataA.numRows);

        DSubmatrixD1 subB = new DSubmatrixD1(B);

        W.col0 = W.row0 = 0;
        Y.row1 = W.row1 = dataA.numRows;
        WTA.row0 = WTA.col0 = 0;

        int start = minDimen - minDimen % blockLength;
        if( start == minDimen )
            start -= blockLength;
        if( start < 0 )
            start = 0;

        // (Q1^T * (Q2^T * (Q3^t * A)))
        for( int i = start; i >= 0; i -= blockLength ) {

            Y.col0 = i;
            Y.col1 = Math.min(Y.col0+blockLength,dataA.numCols);
            Y.row0 = i;
            if( isIdentity )
                subB.col0 = i;
            subB.row0 = i;

            setW();
            WTA.row1 = Y.col1-Y.col0;
            WTA.col1 = subB.col1-subB.col0;
            WTA.original.reshape(WTA.row1,WTA.col1,false);

            // Compute W matrix from reflectors stored in Y
            if( !saveW )
                BlockHouseHolder_DDRB.computeW_Column(blockLength,Y,W,temp, gammas,Y.col0);

            // Apply the Qi to Q
            BlockHouseHolder_DDRB.multTransA_vecCol(blockLength,Y,subB,WTA);
            MatrixMult_DDRB.multPlus(blockLength,W,WTA,subB);
        }
    }
1160214
public static void setIdentity( DMatrixRBlock A )
    {
        int minLength = Math.min(A.numRows,A.numCols);

        CommonOps_DDRM.fill(A, 0);

        int blockLength = A.blockLength;

        for( int i = 0; i < minLength; i += blockLength ) {
            int h = Math.min(blockLength,A.numRows-i);
            int w = Math.min(blockLength,A.numCols-i);

            int index = i*A.numCols + h*i;

            int m = Math.min(h,w);
            for( int k = 0; k < m; k++ ) {
                A.data[index + k*w + k ] = 1;
            }
        }
    }
1075202
public static double findMaxCol(final int blockLength , final DSubmatrixD1 Y , final int col )
    {
        final int width = Math.min(blockLength,Y.col1-Y.col0);

        final double dataY[] = Y.original.data;

        double max = 0;

        for( int i = Y.row0; i < Y.row1; i += blockLength ) {
            int height = Math.min( blockLength , Y.row1 - i );

            int index = i*Y.original.numCols + height*Y.col0 + col;

            if( i == Y.row0 ) {
                index += width*col;
                for( int k = col; k < height; k++ , index += width ) {
                    double v = Math.abs(dataY[index]);
                    if( v > max ) {
                        max = v;
                    }
                }
            } else {
                for( int k = 0; k < height; k++ , index += width ) {
                    double v = Math.abs(dataY[index]);
                    if( v > max ) {
                        max = v;
                    }
                }
            }
        }

        return max;
    }
452242
@Override
    public void invert(DMatrixRMaj A_inv) {
        blockB.reshape(A_inv.numRows,A_inv.numCols,false);

        alg.invert(blockB);

        MatrixOps_DDRB.convert(blockB,A_inv);
    }
964415
@Override
    public boolean setA(DMatrixRMaj A) {
        blockA.reshape(A.numRows,A.numCols,false);
        MatrixOps_DDRB.convert(A,blockA);

        return alg.setA(blockA);
    }
1316486
public void decompose( SimpleMatrix A ) {

        this.QR = A.copy();

        int N = Math.min(A.numCols(),A.numRows());
        gammas = new double[ A.numCols() ];

        for( int i = 0; i < N; i++ ) {
            // use extract matrix to get the column that is to be zeroed
            SimpleMatrix v = QR.extractMatrix(i, END,i,i+1);
            double max = v.elementMaxAbs();

            if( max > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                v = v.divide(max);

                // compute the magnitude of the vector
                double tau = v.normF();

                if( v.get(0) < 0 )
                    tau *= -1.0;

                double u_0 = v.get(0) + tau;
                double gamma = u_0/tau;

                v = v.divide(u_0);
                v.set(0,1.0);

                // extract the submatrix of A which is being operated on
                SimpleMatrix A_small = QR.extractMatrix(i,END,i,END);

                // A = (I - &gamma;*u*u<sup>T</sup>)A
                A_small = A_small.plus(-gamma,v.mult(v.transpose()).mult(A_small));

                // save the results
                QR.insertIntoThis(i,i,A_small);
                QR.insertIntoThis(i+1,i,v.extractMatrix(1,END,0,1));

                // Alternatively, the two lines above can be replaced with in-place equations
                // READ THE JAVADOC TO UNDERSTAND HOW THIS WORKS!
//                QR.equation("QR(i:,i:) = A","QR",i,"i",A_small,"A");
//                QR.equation("QR((i+1):,i) = v(1:,0)","QR",i,"i",v,"v");

                // save gamma for recomputing Q later on
                gammas[i] = gamma;
            }
        }
    }
303892
public static void constraintMatrix3x3a( DMatrixRMaj L_3x6 , DMatrixRMaj L_3x3 ) {

		int index = 0;
		for( int i = 0; i < 3; i++ ) {
			L_3x3.data[index++] = L_3x6.get(i,0);
			L_3x3.data[index++] = L_3x6.get(i,1);
			L_3x3.data[index++] = L_3x6.get(i,2);
		}
	}
435553
public void gaussianElimination(int[][] matrix) {
    assert(matrix != null && matrix.length > 0 && matrix[0].length > 0
           && matrix.length < matrix[0].length);
    int height = matrix.length;
    int width = matrix[0].length;
    for (int i = 0; i < height; i++) {
      boolean pivotFound = false;
      // scan the column for a nonzero pivot and swap it to the diagonal
      for (int j = i; j < height; j++) {
        if (matrix[i][j] != 0) {
          int[] tmp = matrix[i];
          matrix[i] = matrix[j];
          matrix[j] = tmp;
          pivotFound = true;
          break;
        }
      }
      if (!pivotFound) {
        continue;
      }
      int pivot = matrix[i][i];
      for (int j = i; j < width; j++) {
        matrix[i][j] = divide(matrix[i][j], pivot);
      }
      for (int j = i + 1; j < height; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
    for (int i = height - 1; i >=0; i--) {
      for (int j = 0; j < i; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
  }
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;
            }
        }
    }
533367
public static void mult( double a[] , int offsetA ,
                             DMatrixSparseCSC B ,
                             double c[] , int offsetC )
    {
        if( a.length-offsetA < B.numRows)
            throw new IllegalArgumentException("Length of 'a' isn't long enough");
        if( c.length-offsetC < B.numCols)
            throw new IllegalArgumentException("Length of 'c' isn't long enough");

        for (int k = 0; k < B.numCols; k++) {
            int idx0 = B.col_idx[k  ];
            int idx1 = B.col_idx[k+1];

            double sum = 0;
            for (int indexB = idx0; indexB < idx1; indexB++) {
                sum += a[offsetA+B.nz_rows[indexB]]*B.nz_values[indexB];
            }
            c[offsetC+k] = sum;
        }
    }
1367612
public static void multAddTransAB(DMatrixSparseCSC A , DMatrixRMaj B , DMatrixRMaj C )
    {
        if( A.numCols != C.numRows || B.numRows != C.numCols )
            throw new IllegalArgumentException("Inconsistent matrix shapes. "+stringShapes(A,B,C));

        ImplSparseSparseMult_DSCC.multAddTransAB(A,B,C);
    }
780464
@Override
    public boolean setA(DMatrixRBlock A) {
        if( A.numRows < A.numCols )
            throw new IllegalArgumentException("Number of rows must be more than or equal to the number of columns.  " +
                    "Can't solve an underdetermined system.");

        if( !decomposer.decompose(A))
            return false;

        this.QR = decomposer.getQR();

        return true;
    }