sdpa_jordan.cpp

/* -------------------------------------------------------------

This file is a component of SDPA
Copyright (C) 2004 SDPA Project

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA

------------------------------------------------------------- */

#include <sdpa_jordan.h>
#include <sdpa_dataset.h>

namespace sdpa {

mpf_class Jal::trace(DenseLinearSpace& aMat)
{
  mpf_class ret = 0.0;

  for (int l=0; l<aMat.SDP_nBlock; ++l) {
    mpf_class* target = aMat.SDP_block[l].de_ele;
    int size = aMat.SDP_block[l].nRow;
    for (int j=0; j<size; ++j) {
	  ret += target[j*size+j];
	}
  }

  for (int l=0; l<aMat.SOCP_nBlock; ++l) {
    rError("dataset:: current version do not support SOCP");
  }

  for (int l=0; l<aMat.LP_nBlock; ++l) {
    ret += aMat.LP_block[l];
  }

  return ret;
}

// calculate the minimum eigen value of lMat*xMat*(lMat^T)
// lMat is lower triangular¡¢xMat is symmetric
// block size > 20   : Lanczos method
// block size <= 20  : QR method
// QR method: workVec is temporary space and needs
//            3*xMat.nRow-1 length memory.
mpf_class Jal::getMinEigen(DenseLinearSpace& lMat,
		       DenseLinearSpace& xMat,
		       WorkVariables& work)
{
  mpf_class min = 1.0E50;
  mpf_class value;

  // for SDP
  for (int l=0; l<xMat.SDP_nBlock; ++l) {
    if (xMat.SDP_block[l].nRow > 20){ // use Lanczos method
      value = Lal::getMinEigen(lMat.SDP_block[l],
			       xMat.SDP_block[l],
			       work.DLS1.SDP_block[l],
			       work.SDP_BV1.ele[l],
			       work.SDP_BV2.ele[l],
			       work.SDP_BV3.ele[l],
			       work.SDP_BV4.ele[l],
			       work.SDP_BV5.ele[l],
			       work.SDP_BV6.ele[l],
			       work.SDP_BV7.ele[l],
			       work.SDP_BV8.ele[l],
			       work.SDP_BV9.ele[l],
			       work.SDP2_BV1.ele[l]);
    } else { // use QR method
      Lal::let(work.DLS2.SDP_block[l],'=',xMat.SDP_block[l],'T',lMat.SDP_block[l]);
      Lal::let(work.DLS1.SDP_block[l],'=',lMat.SDP_block[l],'*',work.DLS2.SDP_block[l]);
      Lal::getMinEigenValue(work.DLS1.SDP_block[l],work.SDP_BV1.ele[l],work.SDP2_BV1.ele[l]);
      value = work.SDP_BV1.ele[l].ele[0];
    }
    if (value < min) {
      min = value;
    }
  } // end of 'for (int l)'

  // for SOCP
  for (int l=0; l<xMat.SOCP_nBlock; ++l) {
    rError("getMinEigen:: current version does not support SOCP");
  }

// for  LP
  for (int l=0; l<xMat.LP_nBlock; ++l) {
    value = xMat.LP_block[l] * lMat.LP_block[l] * lMat.LP_block[l];
    if (value < min) {
      min = value;
    }
  } // end of 'for (int l)'

  return min;
};


  // calculate the minimum eigen value of xMat by QR method.
mpf_class Jal::getMinEigen(DenseLinearSpace& xMat,
                        WorkVariables& work)
{
  mpf_class min = 1.0E50;
  mpf_class value;

  work.DLS1.copyFrom(xMat);

  // for SDP
  for (int l=0; l<xMat.SDP_nBlock; ++l) {
    Lal::getMinEigenValue(work.DLS1.SDP_block[l],work.SDP_BV1.ele[l],work.SDP2_BV1.ele[l]);
    value = work.SDP_BV1.ele[l].ele[0];

    if (value < min) {
      min = value;
    }
  } // end of 'for (int l)'

  // for SOCP
  for (int l=0; l<xMat.SOCP_nBlock; ++l) {
    rError("getMinEigen:: current version does not support SOCP");
  }

// for  LP
  for (int l=0; l<xMat.LP_nBlock; ++l) {
    value = xMat.LP_block[l];
    if (value < min) {
      min = value;
    }
  } // end of 'for (int l)'

  return min;
};


bool Jal::getInvChol(DenseLinearSpace& invCholMat,
		     DenseLinearSpace& aMat,
		     DenseLinearSpace& workMat)
{
  // for SDP
  if (workMat.SDP_nBlock!=aMat.SDP_nBlock
      || invCholMat.SDP_nBlock!=aMat.SDP_nBlock) {
    rError("getInvChol:: different memory size");
  }
  for (int l=0; l<aMat.SDP_nBlock; ++l) {
    if (Lal::getCholesky(workMat.SDP_block[l],aMat.SDP_block[l]) == false) {
      return false;
    }
    Lal::getInvLowTriangularMatrix(invCholMat.SDP_block[l],
				   workMat.SDP_block[l]);
  }

  // for SOCP
  for (int l=0; l<aMat.SOCP_nBlock; ++l) {
    rError("no support for SOCP");
  }

  // fo LP
  if (invCholMat.LP_nBlock!=aMat.LP_nBlock) {
    rError("getInvChol:: different memory size");
  }
  for (int l=0; l<aMat.LP_nBlock; ++l) {
    if (aMat.LP_block[l] < 0.0){
      return false;
    }
    invCholMat.LP_block[l] = 1.0 / sqrt(aMat.LP_block[l]);
  }

  return _SUCCESS;
}

bool Jal::getInvCholAndInv(DenseLinearSpace& invCholMat,
			   DenseLinearSpace& inverseMat,
			   DenseLinearSpace& aMat,
			   DenseLinearSpace& workMat)
{
  mpf_class MONE = 1.0;
  if (getInvChol(invCholMat, aMat, workMat)== false) {
    return FAILURE;
  }

  for (int l=0; l<aMat.SDP_nBlock; ++l) {
    inverseMat.SDP_block[l].copyFrom(invCholMat.SDP_block[l]);
    Rtrmm ("Left","Lower","Transpose","NonUnitDiag",
	   invCholMat.SDP_block[l].nRow,
	   invCholMat.SDP_block[l].nCol,
	   MONE,
	   invCholMat.SDP_block[l].de_ele,
	   invCholMat.SDP_block[l].nRow,
	   inverseMat.SDP_block[l].de_ele,
           inverseMat.SDP_block[l].nRow);
  }
  for (int l=0; l<aMat.SOCP_nBlock; ++l) {
    rError("rNewton:: we don't make this ruoutin");
    
  }
  for (int l=0; l<aMat.LP_nBlock; ++l) {
    inverseMat.LP_block[l] = 1.0 / aMat.LP_block[l];
  }
  return _SUCCESS;
}


bool Jal::multiply(DenseLinearSpace& retMat,
		   DenseLinearSpace& aMat,
		   DenseLinearSpace& bMat,
		   mpf_class* scalar)
{
  bool total_judge = _SUCCESS;

  // for SDP
  if (retMat.SDP_nBlock!=aMat.SDP_nBlock 
      || retMat.SDP_nBlock!=bMat.SDP_nBlock) {
    rError("multiply:: different nBlock size");
  }
  for (int l=0; l<retMat.SDP_nBlock; ++l) {
    bool judge = Lal::multiply(retMat.SDP_block[l],aMat.SDP_block[l],
			       bMat.SDP_block[l],scalar);
    if (judge == FAILURE) {
      total_judge = FAILURE;
    }
  }

  // for SOCP
#if 0
  if (retMat.SOCP_nBlock!=aMat.SOCP_nBlock 
      || retMat.SOCP_nBlock!=bMat.SOCP_nBlock) {
    rError("multiply:: different nBlock size");
  }
  for (int l=0; l<retMat.SOCP_nBlock; ++l) {
    bool judge = Lal::multiply(retMat.SOCP_block[l],aMat.SOCP_block[l],
			       bMat.SOCP_block[l],scalar);
    if (judge == FAILURE) {
      total_judge = FAILURE;
    }
  }
#endif

  // for LP
  if (retMat.LP_nBlock!=aMat.LP_nBlock 
      || retMat.LP_nBlock!=bMat.LP_nBlock) {
    rError("multiply:: different nBlock size");
  }
  for (int l=0; l<retMat.LP_nBlock; ++l) {
    if (scalar == NULL) {
      retMat.LP_block[l] = aMat.LP_block[l] * bMat.LP_block[l];
    } else {
      retMat.LP_block[l] = aMat.LP_block[l] * bMat.LP_block[l] * (*scalar);
    }
  }

  return total_judge;
}

#if 0
// CAUTION!!! We don't initialize retMat to zero matrix for efficiently.
bool Jal::multiply(DenseLinearSpace& retMat,
		   SparseLinearSpace& aMat,
		   DenseLinearSpace& bMat,
		   mpf_class* scalar)
{
  bool total_judge = _SUCCESS;

  // for SDP
  for (int l=0; l<aMat.SDP_sp_nBlock; ++l) {
    int index = aMat.SDP_sp_index[l];
    bool judge = Lal::multiply(retMat.SDP_block[index],aMat.SDP_sp_block[l],
			       bMat.SDP_block[index],scalar);
    if (judge == FAILURE) {
      total_judge = FAILURE;
    }
  }

  // for SOCP
  for (int l=0; l<aMat.SOCP_sp_nBlock; ++l) {
    int index = aMat.SOCP_sp_index[l];
    bool judge = Lal::multiply(retMat.SOCP_block[index],aMat.SOCP_sp_block[l],
			       bMat.SOCP_block[index],scalar);
    if (judge == FAILURE) {
      total_judge = FAILURE;
    }
  }

  // for LP
  for (int l=0; l<aMat.LP_sp_nBlock; ++l) {
    int index = aMat.LP_sp_index[l];
    if (scalar == NULL) {
      retMat.LP_block[index] = 
	aMat.LP_sp_block[l] *  bMat.LP_block[index];
    } else {
      retMat.LP_block[index] = 
	aMat.LP_sp_block[l] *  bMat.LP_block[index] * (*scalar);
    }
  }

  return total_judge;
}

// CAUTION!!! We don't initialize retMat to zero matrix for efficiently.
bool Jal::multiply(DenseLinearSpace& retMat,
		   DenseLinearSpace& aMat,
		   SparseLinearSpace& bMat,
		   mpf_class* scalar )
{
  bool total_judge = _SUCCESS;

  // for SDP
  for (int l=0; l<bMat.SDP_sp_nBlock; ++l) {
    int index = bMat.SDP_sp_index[l];
    bool judge = Lal::multiply(retMat.SDP_block[index],aMat.SDP_block[index],
			       bMat.SDP_sp_block[l],scalar);
    if (judge == FAILURE) {
      total_judge = FAILURE;
    }
  }

  // for SOCP
#if 0
  for (int l=0; l<bMat.SOCP_sp_nBlock; ++l) {
    int index = bMat.SOCP_sp_index[l];
    bool judge = Lal::multiply(retMat.SOCP_block[index],aMat.SOCP_block[index],
			       bMat.SOCP_sp_block[l],scalar);
    if (judge == FAILURE) {
      total_judge = FAILURE;
    }
  }
#endif

  // for LP
  for (int l=0; l<bMat.LP_sp_nBlock; ++l) {
    int index = bMat.LP_sp_index[l];
    if (scalar == NULL) {
      retMat.LP_block[index] = 
	aMat.LP_block[index] * bMat.LP_sp_block[l];
    } else {
      retMat.LP_block[index] = 
	aMat.LP_block[index] * bMat.LP_sp_block[l] * (*scalar);
    }
  }

  return total_judge;
}
#endif  

  //  retMat = (A * B + B * A)/2
bool Jal::jordan_product(DenseLinearSpace& retMat,
			 DenseLinearSpace& aMat,
			 DenseLinearSpace& bMat)
{
  multiply(retMat,aMat,bMat);
  Lal::getSymmetrize(retMat);
  return _SUCCESS;
}

  //  retMat = A * B
bool Jal::ns_jordan_product(DenseLinearSpace& retMat,
			    DenseLinearSpace& aMat,
			    DenseLinearSpace& bMat)
{
  multiply(retMat,aMat,bMat);
  return _SUCCESS;
}
  
  //  retMat = A * B * A 
bool Jal::jordan_quadratic_product(DenseLinearSpace& retMat,
				   DenseLinearSpace& aMat,
				   DenseLinearSpace& bMat,
				   DenseLinearSpace& workMat)
{
  multiply(workMat,aMat,bMat);
  multiply(retMat,workMat,aMat);
  return _SUCCESS;
}
  
  //  retMat = (A * B * C + C * B * A)/2
bool Jal::jordan_triple_product(DenseLinearSpace& retMat,
				DenseLinearSpace& aMat,
				DenseLinearSpace& bMat,
				DenseLinearSpace& cMat,
				DenseLinearSpace& workMat)
{
  multiply(workMat,aMat,bMat);
  multiply(retMat,workMat,cMat);
  Lal::getSymmetrize(retMat);
  return _SUCCESS;

}

//  retMat = A * B * C
bool Jal::ns_jordan_triple_product(DenseLinearSpace& retMat,
				   DenseLinearSpace& aMat,
				   DenseLinearSpace& bMat,
				   DenseLinearSpace& cMat,
				   DenseLinearSpace& workMat)
{
  multiply(workMat,aMat,bMat);
  multiply(retMat,workMat,cMat);
  return _SUCCESS;
}

}