private void removeDuplicateRow(DoubleMatrix1D c,
			DoubleMatrix2D A, DoubleMatrix1D b, 
			DoubleMatrix1D lb, DoubleMatrix1D ub, 
			DoubleMatrix1D ylb, DoubleMatrix1D yub,
			DoubleMatrix1D zlb, DoubleMatrix1D zub) {
		//the position 0 is for empty rows, 1 is for row singleton and 2 for row doubleton
		int startingLength = 3;
		for(int i=startingLength; i<vRowLengthMap.length; i++){
			int[] vRowLengthMapI = vRowLengthMap[i];
			if(vRowLengthMapI == null || vRowLengthMapI.length < 1){
				//no rows has this number of nz
				continue;
			}
			
			boolean stop = false;
			for(int j=0; !stop && j<vRowLengthMapI.length; j++){
				short prow = (short)vRowLengthMapI[j];//the row of A that has this number of nz
				if(vRowPositions[prow].length==0 || prow < nOfSlackVariables){
					//the upper left part of A is diagonal if there are the slack variables:
					//there is no sparsity superset possible 
					continue;
				}
				short[] vRowPositionsProw = vRowPositions[prow];
				if(vRowPositionsProw.length != i){
					log.debug("Row "+prow+" has an unexpected number of nz: expected " + i + " but is " + vRowPositionsProw.length);
					throw new IllegalStateException();
				}
				for(int si=i; !stop && si<vRowLengthMap.length; si++){
					//look into rows with superset sparsity pattern
					int[] vRowLengthMapSI = vRowLengthMap[si];
					if(vRowLengthMapSI == null || vRowLengthMapSI.length < 1){
						continue;
					}
					for(int sj=0; sj<vRowLengthMapSI.length; sj++){
						if(si==i && sj <= j){
							continue;//look forward, not behind
						}
						short srow = (short)vRowLengthMapSI[sj];
						if(vRowPositions[srow].length==0){
							continue;
						}
						short[] vRowPositionsSrow = vRowPositions[srow];
						//same sparsity pattern?
						if(isSubsetSparsityPattern(vRowPositionsProw, vRowPositionsSrow)){
							log.debug("found superset sparsity pattern: row " + prow + " contained in row " + srow);
							
							//look for the higher number of coefficients that can be deleted
							Map<Double, List<Integer>> coeffRatiosMap = new HashMap<Double, List<Integer>>();
							for(short k=0; k<vRowPositionsProw.length; k++){
								short col = vRowPositionsProw[k];
								double APRL = A.getQuick(prow, col); 
								double ASRL = A.getQuick(srow, col);
								double ratio = -ASRL/APRL;
								//put the ratio and the column index in the map
								boolean added = false;
								for(Double keyRatio : coeffRatiosMap.keySet()){
									if(isZero(ratio - keyRatio)){
										coeffRatiosMap.get(keyRatio).add((int)col);
										added = true;
										break;
									}
								}
								if(!added){
									List<Integer> newList = new ArrayList<Integer>();
									newList.add((int)col);
									coeffRatiosMap.put(ratio, newList);
								}
							}
							//take the ratio(s) with the higher number of column indexes
							int maxNumberOfColumn = -1;
							List<Integer> candidatedColumns = null;
							for(Double keyRatio : coeffRatiosMap.keySet()){
								int size = coeffRatiosMap.get(keyRatio).size(); 
								if(size > maxNumberOfColumn){
									maxNumberOfColumn = size;
									candidatedColumns = coeffRatiosMap.get(keyRatio);
								}else if(size == maxNumberOfColumn){
									candidatedColumns.addAll(coeffRatiosMap.get(keyRatio));
								}
							}

							//look for the position with less column fill in
							short lessFilledColumn = -1;//cannot be greater
							int lessFilledColumnLength = this.originalMeq + 1;//cannot be greater
							for(short k=0; k<candidatedColumns.size(); k++){
								short col = candidatedColumns.get(k).shortValue();
								if(vColPositions[col].length>1 && vColPositions[col].length<lessFilledColumnLength ){
									lessFilledColumn = col;
									lessFilledColumnLength = vColPositions[col].length;
								}
							}
							log.debug("less filled column (" + lessFilledColumn +"): length=" + lessFilledColumnLength);
							double APRL = A.getQuick(prow, lessFilledColumn); 
							double ASRL = A.getQuick(srow, lessFilledColumn);
							double alpha = -ASRL/APRL;
							
							b.setQuick(srow, b.getQuick(srow) + alpha*b.getQuick(prow));
							//substitute A[prow] with A[prow] * alpha*A[row] for every nz entry of A[row]
							for(short t=0; t<vRowPositionsProw.length; t++){
								short cc = vRowPositionsProw[t];
								double nv = 0.;
								if(cc!=lessFilledColumn){
									nv = A.getQuick(srow, cc) + alpha*A.getQuick(prow, cc);
								}
								A.setQuick(srow, cc, nv);
								if(isZero(nv)){
									vRowPositions[srow] = removeElementFromSortedArray(vRowPositions[srow], cc);
									changeColumnsLengthPosition(cc, vColPositions[cc].length, vColPositions[cc].length-1);
									vColPositions[cc] = removeElementFromSortedArray(vColPositions[cc], srow);
									changeRowsLengthPosition(srow, vRowPositions[srow].length+1, vRowPositions[srow].length);
									A.setQuick(srow, cc, 0.);
								}
							}
							this.someReductionDone = true;
							stop = true;
							i=startingLength-1;//restart, ++ comes from the for loop
							break;
						}
					}
				}
			}
		}
	}
private void removeDuplicateRow ( DoubleMatrix1D c , DoubleMatrix2D A , DoubleMatrix1D b , DoubleMatrix1D lb , DoubleMatrix1D ub , DoubleMatrix1D ylb , DoubleMatrix1D yub , DoubleMatrix1D zlb , DoubleMatrix1D zub ) { int startingLength = 3 ; for ( int i = startingLength ; i < vRowLengthMap . length ; i ++ ) { int [ ] vRowLengthMapI = vRowLengthMap [ i ] ; if ( vRowLengthMapI == null || vRowLengthMapI . length < 1 ) { continue ; } boolean stop = false ; for ( int j = 0 ; ! stop && j < vRowLengthMapI . length ; j ++ ) { short prow = ( short ) vRowLengthMapI [ j ] ; if ( vRowPositions [ prow ] . length == 0 || prow < nOfSlackVariables ) { continue ; } short [ ] vRowPositionsProw = vRowPositions [ prow ] ; if ( vRowPositionsProw . length != i ) { log . debug ( "Row " + prow + " has an unexpected number of nz: expected " + i + " but is " + vRowPositionsProw . length ) ; throw new IllegalStateException ( ) ; } for ( int si = i ; ! stop && si < vRowLengthMap . length ; si ++ ) { int [ ] vRowLengthMapSI = vRowLengthMap [ si ] ; if ( vRowLengthMapSI == null || vRowLengthMapSI . length < 1 ) { continue ; } for ( int sj = 0 ; sj < vRowLengthMapSI . length ; sj ++ ) { if ( si == i && sj <= j ) { continue ; } short srow = ( short ) vRowLengthMapSI [ sj ] ; if ( vRowPositions [ srow ] . length == 0 ) { continue ; } short [ ] vRowPositionsSrow = vRowPositions [ srow ] ; if ( isSubsetSparsityPattern ( vRowPositionsProw , vRowPositionsSrow ) ) { log . debug ( "found superset sparsity pattern: row " + prow + " contained in row " + srow ) ; Map < Double , List < Integer > > coeffRatiosMap = new HashMap < Double , List < Integer > > ( ) ; for ( short k = 0 ; k < vRowPositionsProw . length ; k ++ ) { short col = vRowPositionsProw [ k ] ; double APRL = A . getQuick ( prow , col ) ; double ASRL = A . getQuick ( srow , col ) ; double ratio = - ASRL / APRL ; boolean added = false ; for ( Double keyRatio : coeffRatiosMap . keySet ( ) ) { if ( isZero ( ratio - keyRatio ) ) { coeffRatiosMap . get ( keyRatio ) . add ( ( int ) col ) ; added = true ; break ; } } if ( ! added ) { List < Integer > newList = new ArrayList < Integer > ( ) ; newList . add ( ( int ) col ) ; coeffRatiosMap . put ( ratio , newList ) ; } } int maxNumberOfColumn = - 1 ; List < Integer > candidatedColumns = null ; for ( Double keyRatio : coeffRatiosMap . keySet ( ) ) { int size = coeffRatiosMap . get ( keyRatio ) . size ( ) ; if ( size > maxNumberOfColumn ) { maxNumberOfColumn = size ; candidatedColumns = coeffRatiosMap . get ( keyRatio ) ; } else if ( size == maxNumberOfColumn ) { candidatedColumns . addAll ( coeffRatiosMap . get ( keyRatio ) ) ; } } short lessFilledColumn = - 1 ; int lessFilledColumnLength = this . originalMeq + 1 ; for ( short k = 0 ; k < candidatedColumns . size ( ) ; k ++ ) { short col = candidatedColumns . get ( k ) . shortValue ( ) ; if ( vColPositions [ col ] . length > 1 && vColPositions [ col ] . length < lessFilledColumnLength ) { lessFilledColumn = col ; lessFilledColumnLength = vColPositions [ col ] . length ; } } log . debug ( "less filled column (" + lessFilledColumn + "): length=" + lessFilledColumnLength ) ; double APRL = A . getQuick ( prow , lessFilledColumn ) ; double ASRL = A . getQuick ( srow , lessFilledColumn ) ; double alpha = - ASRL / APRL ; b . setQuick ( srow , b . getQuick ( srow ) + alpha * b . getQuick ( prow ) ) ; for ( short t = 0 ; t < vRowPositionsProw . length ; t ++ ) { short cc = vRowPositionsProw [ t ] ; double nv = 0. ; if ( cc != lessFilledColumn ) { nv = A . getQuick ( srow , cc ) + alpha * A . getQuick ( prow , cc ) ; } A . setQuick ( srow , cc , nv ) ; if ( isZero ( nv ) ) { vRowPositions [ srow ] = removeElementFromSortedArray ( vRowPositions [ srow ] , cc ) ; changeColumnsLengthPosition ( cc , vColPositions [ cc ] . length , vColPositions [ cc ] . length - 1 ) ; vColPositions [ cc ] = removeElementFromSortedArray ( vColPositions [ cc ] , srow ) ; changeRowsLengthPosition ( srow , vRowPositions [ srow ] . length + 1 , vRowPositions [ srow ] . length ) ; A . setQuick ( srow , cc , 0. ) ; } } this . someReductionDone = true ; stop = true ; i = startingLength - 1 ; break ; } } } } } }
NB: for the rows of A that contain the slack variables, there cannot be the same sparsity pattern
(A is diagonal in its right-upper part)