SolverIO.java

/*
 * Created on 2011/02/21
 * Copyright (C) 2011 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.sdpj.iocenter;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.FileReader;
import java.io.IOException;
import java.io.PrintStream;
import java.util.ArrayList;

import org.mklab.mpf.ap.MPFloat;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.algorithm.AverageComplementarity;
import org.mklab.sdpj.algorithm.DirectionParameter;
import org.mklab.sdpj.algorithm.Parameter;
import org.mklab.sdpj.algorithm.ParameterType;
import org.mklab.sdpj.algorithm.Phase;
import org.mklab.sdpj.algorithm.RatioInitResCurrentRes;
import org.mklab.sdpj.algorithm.Residuals;
import org.mklab.sdpj.algorithm.Solution;
import org.mklab.sdpj.algorithm.SolveInfo;
import org.mklab.sdpj.algorithm.StepLength;
import org.mklab.sdpj.datastructure.BlockDenseMatrix;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.DenseMatrix;
import org.mklab.sdpj.datastructure.SparseDenseDiagonal;
import org.mklab.sdpj.datastructure.SparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.problem.Problem;
import org.mklab.sdpj.solver.Sdpj;
import org.mklab.sdpj.tool.Tools;


/**
 * @author wu
 * @version $Revision$, 2011/02/21
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
 */
public class SolverIO<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {

  /** */
  private int pIteration;
  /** */
  private int m;

  /** */
  private int[] blockStruct;
  /** */
  private AverageComplementarity<RS,RM,CS,CM> mu;
  /** */
  private RatioInitResCurrentRes<RS,RM,CS,CM> theta;
  /** */
  private SolveInfo<RS,RM,CS,CM> solveInfo;
  /** */
  private StepLength<RS,RM,CS,CM> alpha;
  /** */
  private DirectionParameter<RS,RM,CS,CM> beta;
  /** */
  private Residuals<RS,RM,CS,CM> currentRes;
  /** */
  private PrintStream Save;
  /** */
  private PrintStream Display;

  /** */
  private Phase<RS,RM,CS,CM> phase;
  /** */
  private Solution<RS,RM,CS,CM> currentPt;
  /** */
  private int nDim;
  /** */
  private Vector<RS,RM,CS,CM> b;
  /** */
  private BlockSparseMatrix<RS,RM,CS,CM> C;
  /** */
  private BlockSparseMatrix<RS,RM,CS,CM>[] A;
  /** */
  private Parameter<RS,RM,CS,CM> param;
  /** */
  private String outFile;
  /** */
  private ParameterType parameterType;
  /** */
  private boolean isParameter;
  /** */
  private String paraFile;
  /** */
  private boolean isDataSparse;
  /** */
  private boolean isInitFile;
  /** */
  private String initFile;
  /** */
  private boolean isInitSparse;
  /** */
  private ArrayList<String> matrixData;

  /** */
  private boolean isVsolution;

  /** */
  private boolean isDisplay;

  /** */
  private boolean isVsdpj;

  /** */
  private Vsolution<RS,RM,CS,CM> vsolution;
  
  /** */
  private RS unit;

  /**
   * 新しく生成された<code>IO</code>オブジェクトを初期化します。
   * 
   * @param Save instance of PrintStream for save
   * @param Display instance of PrintStream for display
   * @param unit unit
   */
  public SolverIO(PrintStream Save, PrintStream Display, RS unit) {
    this.unit = unit;
    this.Save = Save;
    this.Display = Display;
    this.A = null;
    this.alpha = null;
    this.b = null;
    this.beta = null;
    this.C = null;
    this.currentPt = null;
    this.currentRes = null;
    this.m = 0;
    this.mu = null;
    this.nDim = 0;
    this.outFile = ""; //$NON-NLS-1$
    this.param = null;
    this.phase = null;
    this.pIteration = 0;
    this.solveInfo = null;
    this.theta = null;

    this.parameterType = null;
    this.isParameter = false;
    this.paraFile = ""; //$NON-NLS-1$
    this.isDataSparse = false;
    this.isInitFile = false;
    this.initFile = ""; //$NON-NLS-1$
    this.isInitSparse = false;
    this.matrixData = null;
  }

  /**
   * @param m m
   */
  public void setM(int m) {
    this.m = m;
  }

  /**
   * @param outFile name of output file
   */
  public void setOutFile(String outFile) {
    this.outFile = outFile;
  }

  /**
   * @param param parameter file
   */
  public void setParam(Parameter<RS,RM,CS,CM> param) {
    this.param = param;
  }

  /**
   * @param A class BlockSparseMatrix information
   */
  public void setA(BlockSparseMatrix<RS,RM,CS,CM>[] A) {
    this.A = A;
  }

  /**
   * @param C class BlockSparseMatrix information
   */
  public void setC(BlockSparseMatrix<RS,RM,CS,CM> C) {
    this.C = C;
  }

  /**
   * @param b is the class Vector information
   */
  public void setB(Vector<RS,RM,CS,CM> b) {
    this.b = b;
  }

  /**
   * @param nDim  size of matrix, which is a block matrix.
   */
  public void setNDim(int nDim) {
    this.nDim = nDim;
  }

  /**
   * @param currentPt current solution
   */
  public void setCurrentPt(Solution<RS,RM,CS,CM> currentPt) {
    this.currentPt = currentPt;
  }

  /**
   * @param phase information of phase
   */
  public void setPhase(Phase<RS,RM,CS,CM> phase) {
    this.phase = phase;
  }

  /**
   * @param currentRes information of current Residuals
   */
  public void setCurrentRes(Residuals<RS,RM,CS,CM> currentRes) {
    this.currentRes = currentRes;
  }

  /**
   * @param beta instance of class DirectionParameter
   */
  public void setBeta(DirectionParameter<RS,RM,CS,CM> beta) {
    this.beta = beta;
  }

  /**
   * @param blockStruct is the block structure vector
   */
  public void setBlockStruct(int[] blockStruct) {
    this.blockStruct = blockStruct;
  }

  /**
   * @param alpha instance of class StepLength
   */
  public void setAlpha(StepLength<RS,RM,CS,CM> alpha) {
    this.alpha = alpha;
  }

  /**
   * @param solveInfo solution information
   */
  public void setSolveInfo(SolveInfo<RS,RM,CS,CM> solveInfo) {
    this.solveInfo = solveInfo;
  }

  /**
   * @param theta instance of RatioInitResCurrentRes
   */
  public void setTheta(RatioInitResCurrentRes<RS,RM,CS,CM> theta) {
    this.theta = theta;
  }

  /**
   * @param mu information of AverageComplementarity
   */
  public void setMu(AverageComplementarity<RS,RM,CS,CM> mu) {
    this.mu = mu;
  }

  /**
   * @param pIteration information of iteration
   */
  public void setPIteration(int pIteration) {
    this.pIteration = pIteration;
  }

  /**
   * @param parameterType set parameter <code>parameterType</code>
   */
  public void setParameterType(ParameterType parameterType) {
    this.parameterType = parameterType;
  }

  /**
   * @param isParameter set parameter <code>isParameter</code>
   */
  public void setIsParameter(boolean isParameter) {
    this.isParameter = isParameter;
  }

  /**
   * @param paraFile set parameter <code>paraFile</code>
   */
  public void setParaFile(String paraFile) {
    this.paraFile = paraFile;
  }

  /**
   * @param isDataSparse set parameter <code>isDataSparse</code>
   */
  public void setIsDataSparse(boolean isDataSparse) {
    this.isDataSparse = isDataSparse;
  }

  /**
   * @param isInitFile set parameter <code>isInitFile</code>
   */
  public void setIsInitFile(boolean isInitFile) {
    this.isInitFile = isInitFile;
  }

  /**
   * @param initFile set parameter <code>initFile</code>
   */
  public void setInitFile(String initFile) {
    this.initFile = initFile;
  }

  /**
   * @param isInitSparse set parameter <code>isInitSparse</code>
   */
  public void setIsInitSparse(boolean isInitSparse) {
    this.isInitSparse = isInitSparse;
  }

  /**
   * @param matrixData set parameter <code>matrixData</code>
   */
  public void setMatrixData(ArrayList<String> matrixData) {
    this.matrixData = matrixData;

  }

  /**
   * components: pIteration, mu, theta, solveInfo, alpha, beta, currentRes, fpout, Display.
   */
  @SuppressWarnings("boxing")
  public void printOneIteration() {
    if (Tools.REVERSE_PRIMAL_DUAL == 1) {
      if (this.Display != null && this.isDisplay == true) {
        this.Display.printf("%2d %s %s %s %s %s %s %s %s\n", this.pIteration, this.mu.getValue().toString("%4.2e"), this.theta.getDual().toString("%4.2e"), this.theta.getPrimal().toString("%4.2e"), //$NON-NLS-1$//$NON-NLS-2$//$NON-NLS-3$//$NON-NLS-4$
            this.solveInfo.getObjValDual().unaryMinus().toString("%+7.2e"), //$NON-NLS-1$
            this.solveInfo.getObjValPrimal().unaryMinus().toString("%+7.2e"), this.alpha.getDual().toString("%4.2e"), this.alpha.getPrimal().toString("%4.2e"), this.beta.getValue().toString("%4.2e")); //$NON-NLS-1$ //$NON-NLS-2$//$NON-NLS-3$//$NON-NLS-4$
      }
      if (this.Save != null) {
        this.Save.printf("%2d %s %s %s %s %s %s %s %s\n", this.pIteration, this.mu.getValue().toString("%4.2e"), this.theta.getDual().toString("%4.2e"), this.theta.getPrimal().toString("%4.2e"), //$NON-NLS-1$//$NON-NLS-2$//$NON-NLS-3$//$NON-NLS-4$
            this.solveInfo.getObjValDual().unaryMinus().toString("%+7.2e"), //$NON-NLS-1$
            this.solveInfo.getObjValPrimal().unaryMinus().toString("%+7.2e"), this.alpha.getDual().toString("%4.2e"), this.alpha.getPrimal().toString("%4.2e"), this.beta.getValue().toString("%4.2e")); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
      }
    } else {
      //this section is no chance to be used!!
      if (this.Display != null && this.isDisplay == true) {
        this.Display.printf("%2d  %s  %s  %s  %s  %s  %s  %s  %s\n", this.pIteration, this.mu.getValue().toString(), this.theta.getPrimal().toString(), this.theta.getDual().toString(), //$NON-NLS-1$
            this.solveInfo.getObjValPrimal().toString(), this.solveInfo.getObjValDual().toString(), this.alpha.getPrimal().toString(), this.alpha.getDual().toString(),
            this.beta.getValue().toString());
      }
      if (this.Save != null) {
        this.Save.printf("%2d %s  %s  %s  %s  %s  %s  %s  %s\n", this.pIteration, this.mu.getValue().toString(), this.theta.getPrimal().toString(), this.theta.getDual().toString(), //$NON-NLS-1$
            this.solveInfo.getObjValPrimal().toString(), this.solveInfo.getObjValDual().toString(), this.alpha.getPrimal().toString(), this.alpha.getDual().toString(),
            this.beta.getValue().toString());
      }
    }
  }

  /**
   * 機種精度を返します。
   * 
   * @param num 精度(10進桁数)
   * @return 機種精度
   */
  private RS createMachineEpsilon(RS num) {
    final RS u = num.createUnit();
    RS value = u.divide(2);

    while (u.add(value).equals(u) == false) {
      value = value.divide(2);
    }

    return value.multiply(2);
  }

  /**
   * components: pIteration, mu, theta, solveInfo, alpha, beta, currentRes, phase, currentPt, cputime, nDim, b, C, A, com, param, fpout, Display.
   * 
   * @throws FileNotFoundException if failure, FileNotFoundException is thrown
   */
  @SuppressWarnings("boxing")
  public void printLastInfo() throws FileNotFoundException {

    //   printOneIteration(fpout, Display);
    //    MPFloat.setDefaultPrecisionByDecimalDigits(77);
    RS mean = ((this.solveInfo.getObjValPrimal().abs()).add(this.solveInfo.getObjValDual().abs())).divide(2);
    RS PDgap = (this.solveInfo.getObjValPrimal().subtract(this.solveInfo.getObjValDual())).abs();
    // double dominator;
    final RS relativeGap;

    if (mean.equals(0)) {
      relativeGap = PDgap;
    } else {
      relativeGap = PDgap.divide(mean);
    }

    final RS gap = this.mu.getValue().multiply(this.nDim);

    // infinity in this case
    final RS digits;
    if (PDgap.isZero()) {
      final RS modifiedPDgap = this.solveInfo.getObjValPrimal().abs().max(1).multiply(createMachineEpsilon(PDgap));
      digits = (modifiedPDgap.divide(mean)).abs().log10().unaryMinus();
    } else {
      digits = (PDgap.divide(mean)).abs().log10().unaryMinus();
    }

    if (this.Display != null && this.isDisplay == true) {
      this.Display.printf("\n"); //$NON-NLS-1$
      this.phase.display(this.Display);
      this.Display.printf("   Iteration = %d\n", this.pIteration); //$NON-NLS-1$
      this.Display.printf("          mu = %s\n", this.mu.getValue().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
      this.Display.printf("relative gap = %s\n", relativeGap.toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
      this.Display.printf("         gap = %s\n", gap.toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
      this.Display.printf("      digits = %s\n", digits.toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$

      if (Tools.REVERSE_PRIMAL_DUAL == 1) {
        this.Display.printf("objValPrimal = %s\n", this.solveInfo.getObjValDual().unaryMinus().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$ 
        this.Display.printf("objValDual   = %s\n", this.solveInfo.getObjValPrimal().unaryMinus().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
        this.Display.printf("p.feas.error = %s\n", this.currentRes.getNormDualMat().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$ 
        this.Display.printf("d.feas.error = %s", this.currentRes.getNormPrimalVec().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$ 
      } else {
        this.solveInfo.display();
        this.currentRes.display();
        this.Display.printf("objValPrimal = %10.16e\n", Double.parseDouble(this.solveInfo.getObjValPrimal().toString())); //$NON-NLS-1$
        this.Display.printf("objValDual   = %10.16e\n", Double.parseDouble(this.solveInfo.getObjValDual().toString())); //$NON-NLS-1$
        this.Display.printf("p.feas.error = %10.16e\n", Double.parseDouble(this.currentRes.getNormPrimalVec().toString())); //$NON-NLS-1$
        this.Display.printf("d.feas.error = %10.16e\n\n", Double.parseDouble(this.currentRes.getNormDualMat().toString())); //$NON-NLS-1$
      }
    }

    if (this.Save != null) {
      this.Save.printf("\n"); //$NON-NLS-1$
      this.phase.display(this.Save);
      this.Save.printf("   Iteration = %d\n", this.pIteration); //$NON-NLS-1$
      this.Save.printf("          mu = %s\n", this.mu.getValue().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
      this.Save.printf("relative gap = %s\n", relativeGap.toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
      this.Save.printf("         gap = %s\n", gap.toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
      this.Save.printf("      digits = %s\n", digits.toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$

      if (Tools.REVERSE_PRIMAL_DUAL == 1) {
        this.Save.printf("objValPrimal = %s\n", this.solveInfo.getObjValDual().unaryMinus().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
        this.Save.printf("objValDual   = %s\n", this.solveInfo.getObjValPrimal().unaryMinus().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
        this.Save.printf("p.feas.error = %s\n", this.currentRes.getNormDualMat().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$
        this.Save.printf("d.feas.error = %s\n\n", this.currentRes.getNormPrimalVec().toString("%4.16e")); //$NON-NLS-1$ //$NON-NLS-2$

      } else {
        this.Save.printf("objValPrimal = %10.10e\n", Double.parseDouble(this.solveInfo.getObjValPrimal().toString())); //$NON-NLS-1$
        this.Save.printf("objValDual   = %10.10e\n", Double.parseDouble(this.solveInfo.getObjValDual().toString())); //$NON-NLS-1$
        this.Save.printf("p.feas.error = %10.10e\n", Double.parseDouble(this.currentRes.getNormPrimalVec().toString())); //$NON-NLS-1$
        this.Save.printf("d.feas.error = %10.10e\n", Double.parseDouble(this.currentRes.getNormDualMat().toString())); //$NON-NLS-1$
      }

      this.Save.printf("\n\nParameters are\n"); //$NON-NLS-1$
      this.param.display(this.Save);

      if (Tools.REVERSE_PRIMAL_DUAL == 1) {
        outputAsMatlabFormat(this.currentPt, this.outFile); //modified by wu,2010/10/15

        this.Save.printf("xVec = \n"); //$NON-NLS-1$
        this.currentPt.getYvec().unaryMinus().display(this.Save);
        this.Save.printf("xMat = \n"); //$NON-NLS-1$
        this.currentPt.getZmat().display(this.Save);

        this.Save.printf("yMat = \n"); //$NON-NLS-1$
        this.currentPt.getXmat().display(this.Save);
      } else {
        this.currentPt.display(this.Save);
      }

      //  final boolean CertificateForResults = true; //control whether to report the certificate of optimality or not.
      if (this.isVsolution) {
        this.runVsolution();
      }
    }
  }

  /**
   * Run class Vsolution to verify the solution of SDPJ
   * 
   * @throws FileNotFoundException If a file error occurs
   */
  private void runVsolution() throws FileNotFoundException {
    Vsolution<RS,RM,CS,CM> vsolution = new Vsolution<>(this.unit);
    vsolution.setDataType(this.unit.createUnit(), MPFloat.getDefaultNumberOfPrecisionDigits());
    if (this.isDataSparse) {
      vsolution.setProblemFileName(this.getFileName(this.outFile) + ".dat-s"); //$NON-NLS-1$
    } else {
      vsolution.setProblemFileName(this.getFileName(this.outFile) + ".dat"); //$NON-NLS-1$
    }
    //  vsdpj.setDataSparse(this.isDataSparse);
    vsolution.setSolutionFileName(this.outFile);
    vsolution.setM(this.m);
    vsolution.setDimension(this.nDim);
    for (int i = 0; i <= this.m; i++) {
      if (i == 0) {
        //F0 must to minus, because exchange primal and dual 
        RS[][] tempC =this.C.getArrayOfElements().clone();
        for (int row = 0; row < tempC.length; row++) {
          for (int col = 0; col < tempC.length; col++) {
            tempC[row][col] = tempC[row][col].multiply(this.unit.createUnit().unaryMinus());
          }
        }
        vsolution.setF(tempC, 0);
      } else {
        vsolution.setF(this.A[i - 1].getArrayOfElements(), i);
      }
    }

    vsolution.setConstantVector(this.b.getElements());
    vsolution.setVectorX(this.currentPt.getYvec().unaryMinus().getElements()); // vector (x) must to minus
    vsolution.setMatX(this.currentPt.getZmat().getArrayOfElements());
    vsolution.setMatY(this.currentPt.getXmat().getArrayOfElements());
    vsolution.setBlockStruct(this.blockStruct);

    if (this.isVsdpj == false) {
      this.Display.println("\n"); //$NON-NLS-1$
      vsolution.setNotShowElapsedTime(true);
      vsolution.run();
      Sdpj.setElapsedTime(vsolution.getElapsedTime());
    } else {
      this.vsolution = vsolution; //save the information of SDP problem
    }

  }

  /**
   * The function <code>getFileName()</code> get a file name.
   * 
   * @param fileName The Name of a file
   * @return a name of the file
   */
  private String getFileName(String fileName) {
    final String target;
    int endIndex = fileName.lastIndexOf("."); //$NON-NLS-1$
    if (endIndex == -1) {
      endIndex = fileName.length();
    }
    target = fileName.substring(0, endIndex);
    return target;
  }

  /**
   * @param currentPt currentPt
   * @param outFile outFile
   * @throws FileNotFoundException in case of fail
   */
  private void outputAsMatlabFormat(Solution<RS,RM,CS,CM> currentPt, String outFile) throws FileNotFoundException {
    final String target;
    int index = outFile.lastIndexOf("/"); //$NON-NLS-1$
    if (index == -1) {
      target = outFile.substring(0);
    } else {
      target = outFile.substring(index + 1);
    }

    String replacement1 = "xVec.txt"; //$NON-NLS-1$
    String outFileName1 = outFile.replaceAll(target, replacement1);
    try (PrintStream fpOutformatlab1 = new PrintStream(new FileOutputStream(new File(outFileName1)))) {
      currentPt.getYvec().unaryMinus().displayAsMatlabformat(fpOutformatlab1);
    }

    String replacement2 = "xMat.txt"; //$NON-NLS-1$
    String outFileName2 = outFile.replaceAll(target, replacement2);
    try (PrintStream fpOutformatlab2 = new PrintStream(new FileOutputStream(new File(outFileName2)))) {
      currentPt.getZmat().displayAsMatlabFormat(fpOutformatlab2);
    }

    String replacement3 = "yMat.txt"; //$NON-NLS-1$
    String outFileName3 = outFile.replaceAll(target, replacement3);
    try (PrintStream fpOutformatlab3 = new PrintStream(new FileOutputStream(new File(outFileName3)))) {
      currentPt.getXmat().displayAsMatlabFormat(fpOutformatlab3);
    }
  }

  /**
   * @param save is a instance of PrintStream
   * @param display is a instance of PrintStream
   */
  public void printHeader(PrintStream save, PrintStream display) {
    if (save != null) {
      save.printf("%2s %5s %10s %8s %8s %9s %9s %8s %7s\n", "n", "mu", "thetaP", "thetaD", "objP", "objD", "alphaP", "alphaD", "beta"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$ //$NON-NLS-7$ //$NON-NLS-8$ //$NON-NLS-9$ //$NON-NLS-10$
    }

    if (display != null && this.isDisplay == true) {
      display.printf("%2s %5s %10s %8s %8s %9s %9s %8s %7s\n", "n", "mu", "thetaP", "thetaD", "objP", "objD", "alphaP", "alphaD", "beta"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$ //$NON-NLS-5$ //$NON-NLS-6$ //$NON-NLS-7$ //$NON-NLS-8$ //$NON-NLS-9$ //$NON-NLS-10$
    }
  }

  /**
   * @param dataFile data file
   * @return problem
   * @throws IOException in case of fail
   */
  public Problem<RS,RM,CS,CM> readFile(String dataFile) throws IOException {
    //RS unit = Tools.getUnitNumber().createUnit();
    Parameter<RS,RM,CS,CM> parameter = new Parameter<>(this.unit);
    parameter.setDefaultParameter(this.parameterType);
    if (this.isParameter) {
      parameter.readFile(new File(this.paraFile));
    }

    try (BufferedReader reader = new BufferedReader(new FileReader(new File(dataFile)))) {

      int m = IO.readInteger(reader, this.Save);
      int nBlock = IO.readInteger(reader.readLine());
      this.blockStruct = IO.readBlockStruct(reader, nBlock);

      int nDim = 0;
      for (int i = 0; i < nBlock; ++i) {
        nDim += Math.abs(this.blockStruct[i]);
      }

      Vector<RS,RM,CS,CM> b = Vector.<RS,RM,CS,CM> read(reader, m, this.unit);

      // count numbers of elements of C and A. Note: C,A must be accessed "twice".
      int[] CNonZeroCount = new int[nBlock];
      int[] ANonZeroCount = new int[nBlock * m];

      // initialize C and A
      readForZeroCount(reader, m, nBlock, this.blockStruct, CNonZeroCount, ANonZeroCount, this.isDataSparse);

      BlockSparseMatrix<RS,RM,CS,CM> C = new BlockSparseMatrix<>(nBlock, this.blockStruct);
      for (int i = 0; i < nBlock; ++i) {
        int size = this.blockStruct[i];
        SparseMatrix<RS,RM,CS,CM> block;
        if (size > 0) {
          block = new SparseMatrix<>(size, size, SparseDenseDiagonal.SPARSE, CNonZeroCount[i], this.unit);
        } else {
          block = new SparseMatrix<>(-size, -size, SparseDenseDiagonal.DIAGONAL, -size, this.unit);
        }
        C.setBlock(i, block);
      }

      BlockSparseMatrix<RS,RM,CS,CM>[] A = new BlockSparseMatrix[m];
      for (int i = 0; i < m; ++i) {
        A[i] = new BlockSparseMatrix<>(nBlock, this.blockStruct);
        for (int j = 0; j < nBlock; ++j) {
          int size = this.blockStruct[j];
          SparseMatrix<RS,RM,CS,CM> block;
          if (size > 0) {
            block = new SparseMatrix<>(size, size, SparseDenseDiagonal.SPARSE, ANonZeroCount[i * nBlock + j], this.unit);
          } else {
            block = new SparseMatrix<>(-size, -size, SparseDenseDiagonal.DIAGONAL, -size, this.unit);
          }
          A[i].setBlock(j, block);
        }
      }

      read(C, A, m, nBlock, this.blockStruct, this.isDataSparse);
      reader.close();

      // if possible , change C and A to Dense

      C.changeToDense();
      for (int i = 0; i < m; ++i) {
        A[i].changeToDense();
      }

      Solution<RS,RM,CS,CM> initPt = new Solution<>(this.unit);

      if (this.isInitFile) {
        initPt.initializeZero(m, nBlock, this.blockStruct);
        if (new File(this.initFile).canRead() == false) {
          Tools.error("Cannot open init file " + this.initFile); //$NON-NLS-1$
        }

        try (BufferedReader fpInit = new BufferedReader(new FileReader(new File(this.initFile)))) {
          read(fpInit, initPt.getXmat(), initPt.getYvec(), initPt.getZmat(), nBlock, this.blockStruct, this.isInitSparse);
          fpInit.close();
        }

        initPt.initializeResetup(nBlock, this.blockStruct);
        this.currentPt = initPt.createClone();
      } else {
        initPt.initialize(m, nBlock, this.blockStruct, parameter.lambdaStar);
        this.currentPt = new Solution<>(m, nBlock, this.blockStruct, parameter.lambdaStar);
      }

      Problem<RS,RM,CS,CM> problem = new Problem<>(m, nBlock, this.blockStruct, A, b, C, parameter, this.currentPt, initPt);
      problem.setNDim(nDim);
      return problem;

    }
  }

  /**
   * isInitDataFile
   * 
   * @param reader is a instance of BufferedReader
   * @param x is a instance of BlockDenseMatrix
   * @param y is a instance of Vector
   * @param z is a instance of BlockDenseMatrix
   * @param nBlock is the number of block
   * @param blockStruct is a integer array
   * @param isSparse is boolean data
   * @throws IOException if failure, IOException is thrown
   */
  private void read(BufferedReader reader, BlockDenseMatrix<RS,RM,CS,CM> x, Vector<RS,RM,CS,CM> y, BlockDenseMatrix<RS,RM,CS,CM> z, int nBlock, int[] blockStruct, boolean isSparse) throws IOException {
    // read initial point
    // y is opposite sign
    //RS unit = Tools.getUnitNumber().createE();
    StringBuffer br1 = new StringBuffer(reader.readLine());
    for (int i = 0; i < y.getDimension(); ++i) {
      RS value = this.unit.valueOf(IO.readOtherStringRealNumber(br1));
      y.setElement(i, value.unaryMinus());
    }

    if (isSparse) {
      // sparse case, z, x in this order

      String line;
      while ((line = reader.readLine()) != null) {
        StringBuffer br2 = new StringBuffer(line);
        int target = IO.readOtherStringInteger(br2);
        int l = IO.readOtherStringInteger(br2);
        int i = IO.readOtherStringInteger(br2);
        int j = IO.readOtherStringInteger(br2);

        RS value = this.unit.valueOf(IO.readOtherStringRealNumber(br2));

        if (blockStruct[l - 1] >= 0) {
          if (target == 1) {
            DenseMatrix<RS,RM,CS,CM> block = z.getBlock(l - 1);
            int nCol = block.getColumnSize();
            block.denseElements[(i - 1) + nCol * (j - 1)] = value;
            block.denseElements[(j - 1) + nCol * (i - 1)] = value;
          } else {
            DenseMatrix<RS,RM,CS,CM> block = x.getBlock(l - 1);
            int nCol = block.getColumnSize();
            block.denseElements[(i - 1) + nCol * (j - 1)] = value;
            block.denseElements[(j - 1) + nCol * (i - 1)] = value;
          }
        } else {
          if (target == 1) {
            DenseMatrix<RS,RM,CS,CM> block = z.getBlock(l - 1);
            block.diagonalElements[i - 1] = value;
          } else {
            DenseMatrix<RS,RM,CS,CM> block = x.getBlock(l - 1);
            block.diagonalElements[i - 1] = value;
          }
        }
      }
    } else {
      // dense case, z, x in this order
      StringBuffer br2 = new StringBuffer(reader.readLine());
      for (int i = 0; i < nBlock; ++i) {
        int size = blockStruct[i];

        DenseMatrix<RS,RM,CS,CM> block = z.getBlock(i);
        if (size > 0) {
          for (int j = 0; j < size * size; ++j) {
            block.denseElements[j] = this.unit.valueOf(IO.readOtherStringRealNumber(br2));
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            block.diagonalElements[j] = this.unit.valueOf(IO.readOtherStringRealNumber(br2));
          }
        }
      }

      StringBuffer br3 = new StringBuffer(reader.readLine());
      for (int i = 0; i < nBlock; ++i) {
        int size = blockStruct[i];
        DenseMatrix<RS,RM,CS,CM> block = x.getBlock(i);
        if (size > 0) {
          for (int j = 0; j < size * size; ++j) {
            block.denseElements[j] = this.unit.valueOf(IO.readOtherStringRealNumber(br3));
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            block.diagonalElements[j] = this.unit.valueOf(IO.readOtherStringRealNumber(br3));
          }
        }
      }
    }
  }

  /**
   * count NonZeroNumber
   * 
   * @param reader is a instance of BufferedReader
   * @param m is the number of matrix
   * @param nBlock is the number of block
   * @param blockStruct is a integer array
   * @param CNonZeroCount is non-zero count of matrix No.0.
   * @param ANonZeroCount is non-zero count of matrix No.x.(x !=0 )
   * @param isSparse is a integer array
   * @throws IOException if failure, IOException is thrown
   */
  private void readForZeroCount(BufferedReader reader, int m, int nBlock, int[] blockStruct, int[] CNonZeroCount, int[] ANonZeroCount, boolean isSparse) throws IOException {
    //RS unit = Tools.getUnitNumber().createE();

    for (int i = 0; i < nBlock; ++i) {
      CNonZeroCount[i] = 0;
    }

    for (int i = 0; i < m; ++i) {
      for (int j = 0; j < nBlock; ++j) {
        ANonZeroCount[i * nBlock + j] = 0;
      }
    }

    if (isSparse) {
      String line;
      while ((line = reader.readLine()) != null) {
        StringBuffer buffer = new StringBuffer(line);
        this.matrixData.add(line);

        int k = IO.readOtherStringInteger(buffer);
        int l = IO.readOtherStringInteger(buffer);
        int i = IO.readOtherStringInteger(buffer);
        int j = IO.readOtherStringInteger(buffer);

        IO.readOtherStringRealNumber(buffer);

        if (k == 0) {
          CNonZeroCount[l - 1]++;
        } else {
          ANonZeroCount[(k - 1) * nBlock + (l - 1)]++;
        }
      }

    } else {
      StringBuffer buffer = new StringBuffer();
      String line;
      while ((line = reader.readLine()) != null) {
        this.matrixData.add(line);
        buffer.append(" " + line); //$NON-NLS-1$
      }

      // k==0
      for (int l = 0; l < nBlock; ++l) {
        int size = blockStruct[l];
        if (size > 0) {
          for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
              RS value = this.unit.valueOf(IO.readOtherStringRealNumber(buffer));
              if (i <= j && value.isZero() == false) {
                CNonZeroCount[l]++;
              }
            }
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            IO.readOtherStringRealNumber(buffer);
            CNonZeroCount[l]++;
          }
        }
      }

      for (int k = 0; k < m; ++k) {
        for (int l = 0; l < nBlock; ++l) {
          int size = blockStruct[l];
          if (size > 0) {
            for (int i = 0; i < size; ++i) {
              for (int j = 0; j < size; ++j) {
                RS value = this.unit.valueOf(IO.readOtherStringRealNumber(buffer));
                if (i <= j && value.isZero() == false) {
                  ANonZeroCount[k * nBlock + l]++;
                }
              }
            }
          } else {
            for (int j = 0; j < -size; ++j) {
              IO.readOtherStringRealNumber(buffer);
              ANonZeroCount[k * nBlock + l]++;
            }
          }
        }
      }
    }
  }

  /**
   * @param C is a instance of BlockSparseMatrix
   * @param A is a instance of BlockSparseMatrix
   * @param m is the number of matrix
   * @param nBlock is the number of block
   * @param argBlockStruct is a integer array
   * @param isSparse is boolean data
   */
  private void read(BlockSparseMatrix<RS,RM,CS,CM> C, BlockSparseMatrix<RS,RM,CS,CM>[] A, int m, int nBlock, int[] argBlockStruct, boolean isSparse) {
    //RS unit = Tools.getUnitNumber().createE();

    for (int i = 0; i < nBlock; ++i) {
      SparseMatrix<RS,RM,CS,CM> block = C.getBlock(i);
      if (block.isSparce()) {
        block.nonZeroCount = 0;
        block.nonZeroEffect = 0;
      }
    }

    for (int i = 0; i < m; ++i) {
      for (int j = 0; j < nBlock; ++j) {
        SparseMatrix<RS,RM,CS,CM> block = A[i].getBlock(j);
        if (block.isSparce()) {
          block.nonZeroCount = 0;
          block.nonZeroEffect = 0;
        }
      }
    }

    if (isSparse) {
      for (int index = 0; index < this.matrixData.size(); index++) {
        StringBuffer buffer = new StringBuffer(this.matrixData.get(index));
        int k = IO.readOtherStringInteger(buffer);
        int l = IO.readOtherStringInteger(buffer);
        int i = IO.readOtherStringInteger(buffer);
        int j = IO.readOtherStringInteger(buffer);

        RS value = this.unit.valueOf(IO.readOtherStringRealNumber(buffer));

        if (k == 0) {
          SparseMatrix<RS,RM,CS,CM> block = C.getBlock(l - 1);
          int count = block.nonZeroCount;
          if (block.isSparce() && count >= block.getNonZeroNumber()) {
            Tools.error("C.ele[" + l + "]" + " is too much element which is assigned"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
          }

          if (C.getTargetBlockStruct(l - 1) > 0) {
            block.rowIndex[count] = i - 1;
            block.columnIndex[count] = j - 1;
            // C must be opposite sign.
            block.sparseElements[count] = value.unaryMinus();
            block.nonZeroCount++;
            if (i == j) {
              block.nonZeroEffect++;
            } else {
              block.nonZeroEffect += 2;
            }
          } else {
            block.diagonalElements[j - 1] = value.unaryMinus();
          }
        } else {
          SparseMatrix<RS,RM,CS,CM> block = A[k - 1].getBlock(l - 1);
          int count = block.nonZeroCount;
          if (block.isSparce() && count >= block.getNonZeroNumber()) {
            Tools.error("A[" + k + "].ele[" + l + "]" + " is too much element which is assigned"); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$ //$NON-NLS-4$
          }

          if (A[k - 1].getTargetBlockStruct(l - 1) > 0) {
            block.rowIndex[count] = i - 1;
            block.columnIndex[count] = j - 1;
            block.sparseElements[count] = value;
            block.nonZeroCount++;
            if (i == j) {
              block.nonZeroEffect++;
            } else {
              block.nonZeroEffect += 2;
            }
          } else {
            block.diagonalElements[j - 1] = value;
          }
        }
      }
    } else {
      StringBuffer buffer = new StringBuffer();
      for (int i = 0; i < this.matrixData.size(); i++) {
        buffer.append(" " + this.matrixData.get(i)); //$NON-NLS-1$
      }

      for (int k = 0; k < nBlock; ++k) {
        SparseMatrix<RS,RM,CS,CM> block = C.getBlock(k);
        int size = argBlockStruct[k];
        if (size > 0) {
          for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
              RS value = this.unit.valueOf(IO.readOtherStringRealNumber(buffer));
              if (i <= j && value.isZero() == false) {
                int count = block.nonZeroCount;
                block.rowIndex[count] = i;
                block.columnIndex[count] = j;
                // C must be opposite sign.
                block.sparseElements[count] = value.unaryMinus();
                block.nonZeroCount++;
                if (i == j) {
                  block.nonZeroEffect++;
                } else {
                  block.nonZeroEffect += 2;
                }
              }
            }
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            // C must be opposite sign.
            block.diagonalElements[j] = this.unit.valueOf(IO.readOtherStringRealNumber(buffer)).unaryMinus();
          }
        }
      }

      for (int k = 0; k < m; ++k) {
        for (int l = 0; l < nBlock; ++l) {
          SparseMatrix<RS,RM,CS,CM> block = A[k].getBlock(l);

          int size = argBlockStruct[l];
          if (size > 0) {
            for (int i = 0; i < size; ++i) {
              for (int j = 0; j < size; ++j) {
                RS value = this.unit.valueOf(IO.readOtherStringRealNumber(buffer));

                if (i <= j && !value.isZero()) {
                  int count = block.nonZeroCount;
                  block.rowIndex[count] = i;
                  block.columnIndex[count] = j;
                  block.sparseElements[count] = value;
                  block.nonZeroCount++;

                  if (i == j) {
                    block.nonZeroEffect++;
                  } else {
                    block.nonZeroEffect += 2;
                  }
                }
              }
            }
          } else {
            for (int j = 0; j < -size; ++j) {
              RS value = this.unit.valueOf(IO.readOtherStringRealNumber(buffer));
              block.diagonalElements[j] = value;
            }
          }
        }
      }
    }

    this.matrixData = new ArrayList<>(1);
  }

  /**
   * Returns phase.
   * 
   * @return phase
   */
  public Phase<RS,RM,CS,CM> getPhase() {
    return this.phase;
  }

  /**
   * vsolutionを設定します。
   * 
   * @param vsolution  vsolution 
   */
  public void setVsolution(boolean vsolution) {
    this.isVsolution = vsolution;
  }

  /**
   * vsolutionを返します。
   * @return vsolution vsolution
   */
  public boolean isVsolution() {
    return this.isVsolution;
  }

  /**
   * isDisplayを設定します。
   * 
   * @param isDisplay isDisplay 
   */
  public void setDisplay(boolean isDisplay) {
    this.isDisplay = isDisplay;
  }

  /**
   * isDisplayを返します。
   * 
   * @return isDisplay isDisplay
   */
  public boolean isDisplay() {
    return this.isDisplay;
  }

  /**
   * @param isVsdpj isVsdpjを設定します。
   */
  public void setVsdpj(boolean isVsdpj) {
    this.isVsdpj = isVsdpj;
  }

  /**
   * isVsdpjを返します。
   * 
   * @return isVsdpj isVsdpj
   */
  public boolean isVsdpj() {
    return this.isVsdpj;
  }

  /**
   * vsolutionを設定します。
   * 
   * @param vsolution vsolution
   */
  public void setVsolution(Vsolution<RS,RM,CS,CM> vsolution) {
    this.vsolution = vsolution;
  }

  /**
   * vsolutionを返します。
   * 
   * @return vsolution vsolution
   */
  public Vsolution<RS,RM,CS,CM> getVsolution() {
    return this.vsolution;
  }

  /**
   * Returns pIteration.
   * 
   * @return pIteration
   */
  public int getpIteration() {
    return this.pIteration;
  }

  /**
   * @return outFile
   */
  public String getOutFile() {
    return this.outFile;
  }
}