text:
#ifndef ALM_I_APG_H
#define ALM_I_APG_H

#include "Matrix.h"
#include "APPROX2.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem of the form f(x)+g(x) under the constraint Mx=c;
// where f(x)=\sum_{j=1}^m lambda_f[j] \phi_j(<A_j,x>)
//and g(x)=sum_{i=1}^n g_i(x_i). We all assume that each \phi_j is 1-smooth.

// Each subproblem solves problem of the form f(x)+<Mx-c, lambda_s>+1/2beta_s\|Mx-c\|^2+g(x) by APG.
// This header file implements ASGARD_DL.

template<typename L, typename D>
class ALM_I_APG: public APPROX2<L, D>
{
private:

  std::vector<D> Au;
  std::vector<D> Az;

  std::vector<D> Mu;
  std::vector<D> Mz;

  std::vector<D> Mx_s;
  std::vector<D> Ax_s;
  std::vector<D> lambda_f;

protected:
  Matrix<L,D> data_A;

  Matrix<L, D> data_M;

  std::vector<D> x_s;

  std::vector<D> lambda_s;

  std::vector<D> M_tlambda_s;

  D beta_s;

  D epsilon_s;

  D m_s;

  D mp_s;

  D m_0;

  L m_1;

  L m_2;

  D c_s;

  D eta;

  D rho;

  D max_Lf_s;

  D max_M_s;

  D Lf;

  D L_m;

  L val_mu_f;
  L val_mu_g;

  std::vector<D> dual_alpha1;   //dual_alpha1[j]=gradient_of_phi_j(<A_j,x>)

  std::vector<D> dual_alpha2;  //dual_alpha2[j]=(<M_j,x>-c_j)

  std::vector<D> baralpha1;  //baralpha1=-sum_{j=1}^m_1 lambda_f[j]*dual_alpha1[j]*A_j

  std::vector<D> baralpha2;   //baralpha2=-sum_{j=1}^m_2 1/beta_s*dual_alpha2[j]*M_j

  D residual1;    //residual1=\dist(-M^\top y, \nabla f(x)+\partial g(x))

  D residual2;    //residual2=\|Mx-c\|

  D function_value;

  L print_every_N_ALM_I_APG;

  D running_time_ALM_I_APG;

  L nb_outer_iters;

  std::vector<D> gradient_of_f;

  ofstream samp_ALM_I_APG;
public:

  virtual inline D gradient_of_phi_j(D, L){return D(NULL);}

  virtual inline D value_of_g_i(D, L){return D(NULL);}
  virtual inline D value_of_phi_j(D, L){return D(NULL);}

  virtual inline D prox_of_g_i(D,D,D, L){return D(NULL);}

  virtual inline D feasible_dual(std::vector<D> &,std::vector<D> &) {return D(NULL);}

  virtual inline D value_of_phistar_i(D,L) {return D(NULL);}
  virtual inline D value_of_g_tilde_star(D, vector<D> &,vector<D> &){return D(NULL);}

  virtual inline void set_matrix_M(){}

  virtual inline void set_matrix_A(){}

  virtual inline D distance_to_subgradient_of_g(){return D(NULL);}
  //virtual inline D distance_to_subgradient_of_g(std::vector<D> &, std::vector<D> & ){return D(NULL);}

  ALM_I_APG(const char* Matrix_file, D val_lambda_f)
  :data_A(), data_M()
  {

  }

  ALM_I_APG()
  :data_A(), data_M()
  {

  }

  L get_nb_features(){
    return data_A.get_d();
  }

   //This function computes \nabla_i f(gamma u +z)
    inline D partial_i_of_f(L i){
      D res=0;
      for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
      {
        L j=data_A.col_idx[k];
        D tmp=lambda_f[j]*gradient_of_phi_j(this->gamma*Au[j]+Az[j], j);
        res+=data_A.A_t[k]*tmp;
      }
      for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
      {
        L j=data_M.col_idx[k];
        D tmp=(this->gamma*Mu[j]+Mz[j]-data_M.b[j])/beta_s+lambda_s[j];
        res+=data_M.A_t[k]*tmp;
      }
      return res;
    }

   //This function computes argmin{x1 t+x2 t^2/2+g_i(x3+t)}
   //inline D compute_prox(D x1, D x2, D x3, L i){
     //return prox_operA_tor(x1,x2,x3,i);
   //}

    void compute_x(){
        for(L i=0;i<this->n;i++){
           this->x[i]=this->gamma*this->u[i]+this->z[i];
           x_s[i]=this->x[i];
          }
        for(L j=0;j<m_2;j++)
            Mx_s[j]=this->gamma*this->Mu[j]+this->Mz[j];
        for(L j=0;j<m_1;j++)
            Ax_s[j]=this->gamma*this->Au[j]+this->Az[j];
    }

   inline void compute_primal_value() {

     D res=0;
     for(L i=0;i<this->n;i++){
       this->x[i]=this->gamma*this->u[i]+this->z[i];
       res+=value_of_g_i(this->x[i],i);
     }

     for(L j=0;j<m_1;j++){
       res+=lambda_f[j]*value_of_phi_j(this->gamma*Au[j]+Az[j],j);
     }
     for(L j=0;j<m_2;j++){
       D tmp=(this->gamma*Mu[j]+Mz[j]-data_M.b[j]+beta_s*lambda_s[j]);
       res+=tmp*tmp/2/beta_s;
     }
     this->primal_value=res;
   }

   inline void compute_dual_value(){
     for(L j=0;j<m_1;j++){
       dual_alpha1[j]=gradient_of_phi_j(this->gamma*Au[j]+Az[j], j);
     }
     for(L j=0;j<m_2;j++){
       dual_alpha2[j]=this->gamma*Mu[j]+Mz[j]-data_M.b[j]+beta_s*lambda_s[j];
     }
     for(L i=0;i<this->n;i++){
       baralpha1[i]=0;
       for(L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++){
         L j=data_A.col_idx[k];
         baralpha1[i]-=dual_alpha1[j]*data_A.A_t[k]*lambda_f[j];
       }
       baralpha2[i]=0;
       for(L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++){
         L j=data_M.col_idx[k];
         baralpha2[i]-=dual_alpha2[j]*data_M.A_t[k]/beta_s;
       }
     }
     D res=0;
     D scal=feasible_dual(baralpha1, baralpha2);
     for(L j=0;j<m_1;j++)
     {
       res-=value_of_phistar_i(scal*dual_alpha1[j],j)*lambda_f[j];
     }
     for(L j=0;j<m_2;j++)
       res-=scal*dual_alpha2[j]*(0.5*scal*dual_alpha2[j]+data_M.b[j]-beta_s*lambda_s[j])/beta_s;

     res-=value_of_g_tilde_star(scal,baralpha1,baralpha2);
     this->dual_value=res;
   }

   inline void compute_function_value() {
     D res=0;
     for(L i=0;i<this->n;i++){
       res+=value_of_g_i(x_s[i],i);
     }
     for(L j=0;j<m_1;j++){
       res+=lambda_f[j]*value_of_phi_j(Ax_s[j],j);
     }
     function_value=res;
   }

   D compute_prox(D x1,D x2,D x3, L i){
       D res=prox_of_g_i(x1+beta_s*(x3-x_s[i]), beta_s+x2,x3,i);
       if(fabs(res)>1e10)
       {cout<<"larger than 1e10:"<<res<<endl;
       cout<<beta_s<<"; "<<x1<<"; "<<x3<<"; "<<x_s[i]<<"; "<<i<<endl;
       D res2=0;
      for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
      {
        L j=data_A.col_idx[k];
        D tmp=lambda_f[j]*gradient_of_phi_j(this->gamma*Au[j]+Az[j], j);
        res2+=data_A.A_t[k]*tmp;
        cout<<"Ax["<<j<<"]="<<this->gamma*Au[j]+Az[j]<<endl;
      }
      for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
      {
        L j=data_M.col_idx[k];
        D tmp=(this->gamma*Mu[j]+Mz[j]-data_M.b[j])/beta_s+lambda_s[j];
        cout<<"Mu["<<j<<"]="<<Mu[j]<<"; gamma="<<this->gamma<<endl;
        res2+=data_M.A_t[k]*tmp;
      }
       }
       return res;
   }

   inline void compute_gradient_norm(){
     std::vector<D> Tx(this->n,0);
     this->do_single_step_prox(Tx);
     D res=0;
     for(L i=0;i<this->n;i++)
      res+=(Tx[i]-this->x[i])*(Tx[i]-this->x[i]);
     this->gradient_norm=sqrt(res);
   }

   inline void set_v()
   {
     this->v.resize(this->n,0);
     D maxv=0;
     D minv=std::numeric_limits<double>::max();
     D sumv=0;
     D sumvi1=0;
     L sumw=0;
     L maxw=0;
     L minw=this->n;
     this->sumofLi=0;
     for(L j=0;j<m_1;j++)
     {
       sumw+=data_A.w_t[j];
       maxw=max(maxw,data_A.w_t[j]);
       minw=min(minw,data_A.w_t[j]);
     }
     cout<<"sumw="<<sumw<<";  maxw="<<maxw<<"; minw="<<minw<<endl;
     for(L i=0;i<this->n;i++)
     {
       D vi=0;
       D vi1=0;
       for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
       {
         L j=data_A.col_idx[k];
         vi+=(1.+(data_A.w_t[j]-1.)*(this->tau-1.)/max(this->n-1.,1.))*data_A.A_t[k]*data_A.A_t[k]*lambda_f[j];
         vi1+=data_A.A_t[k]*data_A.A_t[k]*lambda_f[j];
       }
       for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
       {
         L j=data_M.col_idx[k];
         vi+=(1.+(data_M.w_t[j]-1.)*(this->tau-1.)/max(this->n-1.,1.))*data_M.A_t[k]*data_M.A_t[k]/beta_s;
         vi1+=data_M.A_t[k]*data_M.A_t[k]/beta_s;
       }
       this->v[i]=vi;
       sumv+=vi;
       sumvi1+=vi1;
       if(maxv<vi) maxv=vi;
       if(minv>vi) minv=vi;
     }
     if(this->tau==this->n){
          for(L i=0;i<this->n;i++)
              this->v[i]=sumvi1;
     }
     this->sumofLi=sumvi1;
     cout<<"  max of v: "<<maxv<<" ;  min of v: "<<minv<<" ;  sumofv: "<<sumv<<" sumofLi="<<this->sumofLi<<endl;
   }

   inline void compute_m0(D beta0, D val_eta, D val_rho, L val_tau){
       cout<<"beta0="<<beta0<<"val_eta="<<val_eta<<"; val_rho="<<val_rho<<"; val_tau="<<val_tau<<endl;
     m_s=this->n/val_tau*20;

     //D tmp7=1-val_tau/this->n*sqrt(beta0/(max_Lf_s+max_M_s/beta0+beta0));
     //cout<<"tmp7: "<<tmp7<<"; "<<val_tau/(this->n+0.)*sqrt(beta0/(max_Lf_s+max_M_s/beta0+beta0))<<endl;
     //m_0=(2*log(val_rho)+log(val_eta)+log(2))/log(tmp7);

     cout<<" m_s="<<m_s<<endl;
   }

  inline void rescale_Matrix(){
    Lf=0;
    max_Lf_s=0;
    max_M_s=0;
    L_m=0;
    cout<<"this->n="<<this->n<<endl;
    for(L i=0;i<this->n;i++)
    {
      D vi=0;
      for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
      {
        L j=data_A.col_idx[k];
        Lf+=data_A.A_t[k]*data_A.A_t[k]*lambda_f[j];
        vi+=(1.+(data_A.w_t[j]-1.)*(this->tau-1.)/max(this->n-1.,1.))*data_A.A_t[k]*data_A.A_t[k]*lambda_f[j];
      }
      if(vi>max_Lf_s) max_Lf_s=vi;
      D vi2=0;
      for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
      {
        L j=data_M.col_idx[k];
        L_m+=data_M.A_t[k]*data_M.A_t[k];
        vi2+=(1.+(data_M.w_t[j]-1.)*(this->tau-1.)/max(this->n-1.,1.))*data_M.A_t[k]*data_M.A_t[k];
      }
      if(vi2>max_M_s) {max_M_s=vi2;cout<<vi2<<" ";}
    }
    cout<<"max_M_s"<<max_M_s<<endl;
    D scal=sqrt((Lf+1)/L_m);
    cout<<"scal="<<scal<<endl;
    if(scal>1)
    {
      for(L i=0;i<this->n;i++){
        for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
        {
          data_M.A_t[k]*=scal;
        }
      }
      for(L j=0;j<m_2;j++){
        for (L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++)
        {
          data_M.A[k]*=scal;
        }
      }
      for(L j=0;j<m_2;j++)
        data_M.b[j]*=scal;
      max_M_s*=scal*scal;
      L_m*=scal*scal;
    }
    cout<<"max_M_s="<<max_M_s<<"; L_m="<<L_m<<"; max_Lf_s="<<max_Lf_s<<"; L_f="<<Lf<<endl;

  }

   inline void set_p(){
     this->proba_vector.resize(this->n,(0.0+this->tau)/this->n);
     this->max_p=(0.0+this->tau)/this->n;
   }

   void compute_KKT_residual(){
       D res;
       for(L i=0; i<this->n;i++){
           res=0;
           for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
           {
                L j=data_A.col_idx[k];
                D tmp=lambda_f[j]*gradient_of_phi_j(Ax_s[j], j);
                res+=data_A.A_t[k]*tmp;
           }
           gradient_of_f[i]=res;
       }
       residual1=distance_to_subgradient_of_g();
       res=0;
       for(L j=0;j<m_2;j++)
           res+=(Mx_s[j]-data_M.b[j])*(Mx_s[j]-data_M.b[j]);
       residual2=sqrt(res);
   }

   void testing(){
       std::vector<D> nablaf(this->n);
       D res1=0;
       compute_Mty();
       for(L i=0;i<this->n;i++)
       {
          nablaf[i]=0;
          for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
           {
             L j=data_M.col_idx[k];
             D tmp=lambda_s[j];
             nablaf[i]+=data_M.A_t[k]*tmp;
            }
          res1+=(nablaf[i]-M_tlambda_s[i])*(nablaf[i]-M_tlambda_s[i]);
       }
       //samp_ALM2<<" res1 ="<<sqrt(res1)<<endl;
   }

   inline void update_z_coordinate( L i, D dz){
     this->z[i]+=dz;
     for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
     {
       L j=data_A.col_idx[k];
       Az[j]+=dz*data_A.A_t[k];
     }
     for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
     {
       L j=data_M.col_idx[k];
       Mz[j]+=dz*data_M.A_t[k];
       //if(j==2869) cout<<"gamma="<<this->gamma<<"; "<<dz<<"; "<<Mz[j]<<endl;
     }
   }

    inline void update_x_coordinate( L i, D dx){
     this->x[i]+=dx;
     x_s[i]=this->x[i];
     L j;
     for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
     {
       j=data_M.col_idx[k];
       Mx_s[j]+=dx*data_M.A_t[k];
     }
      for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
     {
       j=data_A.col_idx[k];
       Ax_s[j]+=dx*data_A.A_t[k];
     }
   }

   inline void update_u_coordinate( L i, D du){
     this->u[i]+=du;
     for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
     {
       L j=data_A.col_idx[k];
       Au[j]+=du*data_A.A_t[k];
     }
     for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
     {
       L j=data_M.col_idx[k];
       Mu[j]+=du*data_M.A_t[k];
     }
   }

   void compute_Au(){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Au[j]+=this->u[kj]*data_A.A[k];
       }
   }

   void compute_Mu(){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mu[j]+=this->u[kj]*data_M.A[k];
       }
   }

   void compute_Az(){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Az[j]+=this->z[kj]*data_A.A[k];
       }
   }
   void compute_Mz(){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mz[j]+=this->z[kj]*data_M.A[k];
       }
   }

   void compute_Au(vector<D> & x0){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Au[j]+=x0[kj]*data_A.A[k];
       }

   }

   void compute_Az(vector<D> & x0){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Az[j]+=x0[kj]*data_A.A[k];
       }
   }
   void compute_Mu(vector<D> & x0){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mu[j]+=x0[kj]*data_M.A[k];
       }

   }

   void compute_Mz(vector<D> & x0){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mz[j]+=x0[kj]*data_M.A[k];
       }
   }

     void compute_Mx(vector<D> & x0){
     for(L j=0;j<m_2;j++)
       for(L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++){
         L kj=data_M.row_idx[k];
         Mx_s[j]+=x0[kj]*data_M.A[k];
       }
   }

   void compute_Ax(vector<D> & x0){
     for(L j=0;j<m_1;j++)
       for(L k = data_A.ptr[j]; k < data_A.ptr[j + 1];k++){
         L kj=data_A.row_idx[k];
         Ax_s[j]+=x0[kj]*data_A.A[k];
       }
   }

   void compute_cs(){
       c_s=0;
       for(L j=0;j<m_2;j++)
           c_s+=data_M.b[j]*lambda_s[j];
   }

    void compute_Mty(){
        M_tlambda_s.clear();
        M_tlambda_s.resize(this->n,0);
        for(L i=0;i<this->n;i++)
        {
          for(L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++){
            L j=data_M.col_idx[k];
            M_tlambda_s[i]+=lambda_s[j]*data_M.A_t[k];
           }
        }

    }

    void update_y(){
      D res=0;
      D dl=0;
        for(L j=0;j<m_2;j++)
           {
             dl=(Mx_s[j]-data_M.b[j])/beta_s;
             lambda_s[j]+=dl;
             res+=(Mx_s[j]-data_M.b[j])*(Mx_s[j]-data_M.b[j]);
             L i;
             for (L k = data_M.ptr[j]; k < data_M.ptr[j + 1];k++)
             {
               i=data_M.row_idx[k];
               M_tlambda_s[i]+=dl*data_M.A[k];
             }
            }
       //samp_ALM2<<"lambda_{k+1}-lambda_k"<<sqrt(res)<<endl;
    }

   void Initialize(D beta_0, D epsilon_0,  D val_eta, D val_rho,L val_tau, vector<D> & x0,vector<D> & y0, D val_lambda_f){
        cout<<"start initializing"<<endl;
        set_matrix_M();
        set_matrix_A();

        this->tau=val_tau;
        m_1=data_A.get_n();
        m_2=data_M.get_n();
        cout<<"m_1="<<m_1<<endl;
        cout<<"m_2="<<m_2<<endl;
        lambda_f.resize(m_1,val_lambda_f);
        dual_alpha1.resize(m_1,0);
        dual_alpha2.resize(m_2,0);
        this->n=data_A.nfeatures;
        baralpha1.resize(this->n,0);
        baralpha2.resize(this->n,0);
        gradient_of_f.resize(this->n,0);

       rescale_Matrix();

       beta_s=beta_0;
       //if(max_Lf_s>0)
       //  beta_s=min(beta_0,max_M_s/max_Lf_s);

       epsilon_s=epsilon_0;
       compute_m0(beta_s,val_eta,val_rho,val_tau);

       eta=val_eta;
       rho=val_rho;

       x_s.resize(this->n,0);
       lambda_s.resize(m_2,0);
       for(L i=0;i<this->n;i++)
           x_s[i]=x0[i];

       for(L j=0;j<m_2;j++)
           lambda_s[j]=y0[j];

       M_tlambda_s.resize(this->n,0);

       compute_cs();
       compute_Mty();

       Au.clear();
       Au.resize(m_1,0);
       Az.clear();
       Az.resize(m_1,0);
       Ax_s.clear();
       Ax_s.resize(m_1,0);
       Mu.clear();
       Mu.resize(m_2,0);
       Mz.clear();
       Mz.resize(m_2,0);
       Mx_s.clear();
       Mx_s.resize(m_2,0);

       compute_Az(x0);
       compute_Mz(x0);
       compute_Ax(x0);
       compute_Mx(x0);
   }

   void reset_everything(){
       epsilon_s*=rho*rho*eta;
       m_s/=eta;
       beta_s*=eta;

       compute_cs();
       compute_Mty();

       Au.clear();
       Au.resize(m_1,0);
       Az.clear();
       Az.resize(m_1,0);
       Mu.clear();
       Mu.resize(m_2,0);
       Mz.clear();
       Mz.resize(m_2,0);
       compute_Az(x_s);
       compute_Mz(x_s);
   }

    void APPROX_restart(L K, vector<D> & x0, L val_tau, L eval, L p_N,  L max_nb, D eps, string filename){
        std::vector<D> tmp_x(this->n);
        for(L i=0;i<this->n;i++){
          tmp_x[i]=x0[i];
       }
       L k=0;
       while(k<ceil(max_nb/(K+0.0))){
          k++;
          Au.clear();
          Au.resize(m_1,0);
          Az.clear();
          Az.resize(m_1,0);
          Mu.clear();
          Mu.resize(m_2,0);
          Mz.clear();
          Mz.resize(m_2,0);
          compute_Az(tmp_x);
          compute_Mz(tmp_x);
          cout<<"p_N="<<tmp_x[0]<<" "<<tmp_x[1]<<" "<<tmp_x[2]<<" "<<p_N<<endl;
          this->APPROX_MU(tmp_x, val_tau, 0, 0, eval, p_N, K, eps, filename,1);

           for(L i=0;i<this->n;i++){
             tmp_x[i]=this->x[i];
          }
          cout<<"restarting..."<<k<<"; "<<(max_nb/(K+0.0))<<endl;
          if(this->delta<eps) k=ceil(max_nb/(K+0.0))+1;
        }

  }

   inline void compute_and_record_res(){
        if(nb_outer_iters%print_every_N_ALM_I_APG==0){
          compute_KKT_residual();
          compute_function_value();
          cout<<setprecision(9)<<"Iteration: "<<nb_outer_iters<<"; time="<<running_time_ALM_I_APG<<"; residual1="<<residual1<<" residual2="<<residual2<<"; function value="<<function_value<<endl;
          samp_ALM_I_APG<<setprecision(9)<<nb_outer_iters<<" "<<running_time_ALM_I_APG<<" "<<residual1<<" "<<residual2<<" "<<function_value<<" "<<endl;
        }
   }

   void ALM_I_APG_solve_with_APPROX(D beta_0, D epsilon_0,  D eta, D rho,vector<D> & x0,vector<D> & y0, L val_tau, L max_nb_outer, L p_N_1, L p_N_2, D val_lambda_f,string filename1, string filename2, D time){
      Initialize(beta_0, epsilon_0, eta, rho,val_tau, x0, y0, val_lambda_f);

      nb_outer_iters=0;
      string sampname2="results/APPROXMU_"+filename2;
      this->samp.open(sampname2.c_str());
      filename1="results/ALM2_"+filename1;
      samp_ALM_I_APG.open(filename1.c_str());
      running_time_ALM_I_APG=0;
      print_every_N_ALM_I_APG=p_N_1;
      compute_and_record_res();
      D start;
      start = std::clock();
      std::vector<D> record_x_s(this->n);
      std::vector<D> record_x_sp1(this->n);
      cout<<"m0="<<ceil(m_s/this->n*val_tau)<<endl;
      //D K=ceil(2*sqrt((max_Lf_s+max_M_s)/(beta_s+val_mu_g)+1));
      //cout<<"restart period="<<K<<endl;
      cout<<"ms="<<ceil(m_s/this->n*val_tau)<<"; beta_s="<<beta_s<<"; epsilon_s="<<epsilon_s<<endl;
      this->APPROX_MU(x_s, val_tau, val_mu_f, val_mu_g, 3, p_N_2, ceil(m_s/this->n*val_tau), epsilon_s, filename2,1);
      //APPROX_restart(K, x_s, val_tau, 3, p_N_2, ceil((m_s+mp_s)/this->n*val_tau), epsilon_s, filename2);
      compute_x();
      update_y();
      nb_outer_iters++;
      running_time_ALM_I_APG+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
      compute_and_record_res();
      start = std::clock();
      reset_everything();
      running_time_ALM_I_APG+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
      while(nb_outer_iters<max_nb_outer){
         start = std::clock();
         cout<<"ms="<<ceil(m_s/this->n*val_tau)<<"; beta_s="<<beta_s<<"; epsilon_s"<<epsilon_s<<endl;
         //D K=ceil(2*sqrt((max_Lf_s+max_M_s)/(beta_s+val_mu_g)+1));
         //cout<<"restart period="<<K<<endl;
        this->APPROX_MU(x_s, val_tau, val_mu_f, val_mu_g, 3, p_N_2, ceil(m_s/this->n*val_tau), epsilon_s, filename2,1);
        for(L i=0;i<this->n;i++){
          record_x_s[i]=x_s[i];
        }
         //APPROX_restart(K, x_s, val_tau, 3, p_N_2, ceil((m_s+mp_s)/this->n*val_tau), epsilon_s, filename2);
         compute_x();
         update_y();
         nb_outer_iters++;
         running_time_ALM_I_APG+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
         testing();
         compute_and_record_res();
         start = std::clock();
         reset_everything();
         running_time_ALM_I_APG+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
         if (running_time_ALM_I_APG> time){
         	break;
		 }
      }

   }

};

#endif


text:
/* ode-initval2/modnewton1.c
 *
 * Copyright (C) 2008, 2009, 2010 Tuomo Keskitalo
 *
 * 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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* A modified Newton iteration method for solving non-linear
   equations in implicit Runge-Kutta methods.
*/

/* References:

   Ascher, U.M., Petzold, L.R., Computer methods for ordinary
   differential and differential-algebraic equations, SIAM,
   Philadelphia, 1998. ISBN 0898714125, 9780898714128

   Hairer, E., Wanner, G., Solving Ordinary Differential
   Equations II: Stiff and Differential-Algebraic Problems,
   Springer, 1996. ISBN 3540604529, 9783540604525
*/

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

#include "odeiv_util.h"

typedef struct
{
  /* iteration matrix I - h A (*) J */
  gsl_matrix *IhAJ;

  /* permutation for LU-decomposition */
  gsl_permutation *p;

  /* difference vector for kth Newton iteration */
  gsl_vector *dYk;

  /* scaled dYk with desired error level */
  gsl_vector *dScal;

  /* current Runge-Kutta points (absolute values) */
  double *Yk;

  /* function values at points Yk */
  double *fYk;

  /* right hand side of Newton iteration formula */
  gsl_vector *rhs;

  /* stopping criterion value from previous step */
  double eeta_prev;
}
modnewton1_state_t;

static void *
modnewton1_alloc (size_t dim, size_t stage)
{
  modnewton1_state_t *state =
    (modnewton1_state_t *) malloc (sizeof (modnewton1_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for modnewton1_state_t",
                      GSL_ENOMEM);
    }

  state->IhAJ = gsl_matrix_alloc (dim * stage, dim * stage);

  if (state->IhAJ == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for IhAJ", GSL_ENOMEM);
    }

  state->p = gsl_permutation_alloc (dim * stage);

  if (state->p == 0)
    {
      gsl_matrix_free (state->IhAJ);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for p", GSL_ENOMEM);
    }

  state->dYk = gsl_vector_alloc (dim * stage);

  if (state->dYk == 0)
    {
      gsl_permutation_free (state->p);
      gsl_matrix_free (state->IhAJ);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dYk", GSL_ENOMEM);
    }

  state->dScal = gsl_vector_alloc (dim * stage);

  if (state->dScal == 0)
    {
      gsl_vector_free (state->dYk);
      gsl_permutation_free (state->p);
      gsl_matrix_free (state->IhAJ);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dScal", GSL_ENOMEM);
    }

  state->Yk = (double *) malloc (dim * stage * sizeof (double));

  if (state->Yk == 0)
    {
      gsl_vector_free (state->dScal);
      gsl_vector_free (state->dYk);
      gsl_permutation_free (state->p);
      gsl_matrix_free (state->IhAJ);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for Yk", GSL_ENOMEM);
    }

  state->fYk = (double *) malloc (dim * stage * sizeof (double));

  if (state->fYk == 0)
    {
      free (state->Yk);
      gsl_vector_free (state->dScal);
      gsl_vector_free (state->dYk);
      gsl_permutation_free (state->p);
      gsl_matrix_free (state->IhAJ);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for Yk", GSL_ENOMEM);
    }

  state->rhs = gsl_vector_alloc (dim * stage);

  if (state->rhs == 0)
    {
      free (state->fYk);
      free (state->Yk);
      gsl_vector_free (state->dScal);
      gsl_vector_free (state->dYk);
      gsl_permutation_free (state->p);
      gsl_matrix_free (state->IhAJ);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for rhs", GSL_ENOMEM);
    }

  state->eeta_prev = GSL_DBL_MAX;

  return state;
}

static int
modnewton1_init (void *vstate, const gsl_matrix * A,
                 const double h, const gsl_matrix * dfdy,
                 const gsl_odeiv2_system * sys)
{
  /* Initializes the method by forming the iteration matrix IhAJ
     and generating its LU-decomposition
   */

  modnewton1_state_t *state = (modnewton1_state_t *) vstate;

  gsl_matrix *const IhAJ = state->IhAJ;
  gsl_permutation *const p = state->p;

  const size_t dim = sys->dimension;
  const size_t stage = A->size1;

  state->eeta_prev = GSL_DBL_MAX;

  /* Generate IhAJ */

  {
    size_t i, j, k, m;

    for (i = 0; i < dim; i++)
      for (j = 0; j < dim; j++)
        for (k = 0; k < stage; k++)
          for (m = 0; m < stage; m++)
            {
              size_t x = dim * k + i;
              size_t y = dim * m + j;

              if (x != y)
                gsl_matrix_set (IhAJ, x, y,
                                -h * gsl_matrix_get (A, k, m) *
                                gsl_matrix_get (dfdy, i, j));
              else
                gsl_matrix_set (IhAJ, x, y,
                                1.0 - h * gsl_matrix_get (A, k, m) *
                                gsl_matrix_get (dfdy, i, j));
            }
  }

  /* decompose */

  {
    int signum;
    int s = gsl_linalg_LU_decomp (IhAJ, p, &signum);

    if (s != GSL_SUCCESS)
      return s;
  }

  return GSL_SUCCESS;
}

static int
modnewton1_solve (void *vstate, const gsl_matrix * A,
                  const double c[], const double t, const double h,
                  const double y0[], const gsl_odeiv2_system * sys,
                  double YZ[], const double errlev[])
{
  /* Solves the non-linear equation system resulting from implicit
     Runge-Kutta methods by a modified Newton iteration. The
     modified Newton iteration with this problem is stated in the
     form

     IhAJ * dYk = rhs

     in which IhAJ is the iteration matrix: IhAJ = (I - hA (*) J) in
     which (*) is the Kronecker matrix product (size of IhAJ =
     dim*stage, dim*stage), dYk is the Runge-Kutta point (Y)
     difference vector for kth Newton iteration: dYk = Y(k+1) - Y(k),
     and rhs = Y(k) - y0 - h * sum j=1..stage (a_j * f(Y(k)))

     This function solves dYk by LU-decomposition of IhAJ with partial
     pivoting.
   */

  modnewton1_state_t *state = (modnewton1_state_t *) vstate;

  gsl_matrix *const IhAJ = state->IhAJ;
  gsl_permutation *const p = state->p;
  gsl_vector *const dYk = state->dYk;
  double *const Yk = state->Yk;
  double *const fYk = state->fYk;
  gsl_vector *const rhs = state->rhs;
  double *const eeta_prev = &(state->eeta_prev);
  gsl_vector *const dScal = state->dScal;

  const size_t dim = sys->dimension;
  const size_t stage = A->size1;

  int status = GSL_CONTINUE;    /* Convergence status for Newton iteration */

  /* Set starting values for iteration. Use starting values of Y(k) =
     y0. FIXME: Implement better initial guess described in Hairer and
     Wanner.
   */

  {
    size_t j, m;

    gsl_vector_set_zero (dYk);

    for (j = 0; j < stage; j++)
      for (m = 0; m < dim; m++)
        Yk[j * dim + m] = y0[m];
  }

  /* Loop modified Newton iterations */

  {
    size_t iter = 0;
    size_t j, m, n;

    /* Maximum number of Newton iterations. */
    const size_t max_iter = 7;

    double dScal_norm = 0.0;    /* Newton iteration step length */
    double dScal_norm_prev = 0.0;       /* Previous dScal_norm */

    while (status == GSL_CONTINUE && iter < max_iter)
      {
        iter++;

        /* Update Y(k) and evaluate function */

        for (j = 0; j < stage; j++)
          {
            for (m = 0; m < dim; m++)
              {
                Yk[j * dim + m] += gsl_vector_get (dYk, j * dim + m);
              }

            {
              int s = GSL_ODEIV_FN_EVAL (sys, t + c[j] * h, &Yk[j * dim],
                                         &fYk[j * dim]);
              if (s != GSL_SUCCESS)
                {
                  return s;
                }
            }
          }

        /* Calculate rhs  */

        for (j = 0; j < stage; j++)
          for (m = 0; m < dim; m++)
            {
              double sum = 0;

              for (n = 0; n < stage; n++)
                sum += gsl_matrix_get (A, j, n) * fYk[n * dim + m];

              gsl_vector_set (rhs, j * dim + m,
                              -1.0 * Yk[j * dim + m] + y0[m] + h * sum);
            }

        /* Solve dYk */

        {
          int s = gsl_linalg_LU_solve (IhAJ, p, rhs, dYk);

          if (s != GSL_SUCCESS)
            {
              return s;
            }
        }

        /* Evaluate convergence according to method by Hairer and
           Wanner, section IV.8. The iteration step is normalized by
           desired error level before calculation of the norm, to take
           the tolerance on each component of y into account.
         */

        {
          /* relative change in two Newton iteration steps */
          double theta_k;

          /* transformation of theta_k */
          double eeta_k = 0.0;

          for (j = 0; j < stage; j++)
            for (m = 0; m < dim; m++)
              {
                gsl_vector_set (dScal, j * dim + m,
                                gsl_vector_get (dYk, j * dim + m)
                                / errlev[m]);
              }

          dScal_norm_prev = dScal_norm;
          dScal_norm = gsl_blas_dnrm2 (dScal);

          theta_k = dScal_norm / dScal_norm_prev;

          if (iter > 1)
            {
              /* check for increase in step size, which indicates divergence */

              if (theta_k >= 1.0)
                {
                  return GSL_FAILURE;
                }

              eeta_k = theta_k / (1.0 - theta_k);
            }

          else
            {
              eeta_k = pow (GSL_MAX_DBL (*eeta_prev, GSL_DBL_EPSILON), 0.8);
            }

          *eeta_prev = eeta_k;

          if (eeta_k * dScal_norm < 1.0)
            {
              status = GSL_SUCCESS;
            }
        }
      }
  }

  /* No convergence reached within allowed iterations */

  if (status != GSL_SUCCESS)
    {
      return GSL_FAILURE;
    }

  /* Give solution in YZ */

  else
    {
      size_t j, m;

      for (j = 0; j < stage; j++)
        for (m = 0; m < dim; m++)
          YZ[j * dim + m] =
            Yk[j * dim + m] + gsl_vector_get (dYk, j * dim + m);

      return status;
    }
}

static void
modnewton1_free (void *vstate)
{
  modnewton1_state_t *state = (modnewton1_state_t *) vstate;

  gsl_vector_free (state->rhs);
  free (state->fYk);
  free (state->Yk);
  gsl_vector_free (state->dScal);
  gsl_vector_free (state->dYk);
  gsl_permutation_free (state->p);
  gsl_matrix_free (state->IhAJ);
  free (state);
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* multifit/covar.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Brian Gough
 *
 * 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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlin.h>

/* Compute the covariance matrix

   cov = inv (J^T J)

   by QRP^T decomposition of J
*/

int
gsl_multifit_covar (const gsl_matrix * J, double epsrel, gsl_matrix * covar)
{
  double tolr;

  size_t i, j, k;
  size_t kmax = 0;

  gsl_matrix * r;
  gsl_vector * tau;
  gsl_vector * norm;
  gsl_permutation * perm;

  size_t m = J->size1, n = J->size2 ;

  if (m < n)
    {
      GSL_ERROR ("Jacobian be rectangular M x N with M >= N", GSL_EBADLEN);
    }

  if (covar->size1 != covar->size2 || covar->size1 != n)
    {
      GSL_ERROR ("covariance matrix must be square and match second dimension of jacobian", GSL_EBADLEN);
    }

  r = gsl_matrix_alloc (m, n);
  tau = gsl_vector_alloc (n);
  perm = gsl_permutation_alloc (n) ;
  norm = gsl_vector_alloc (n) ;

  {
    int signum = 0;
    gsl_matrix_memcpy (r, J);
    gsl_linalg_QRPT_decomp (r, tau, perm, &signum, norm);
  }

  /* Form the inverse of R in the full upper triangle of R */

  tolr = epsrel * fabs(gsl_matrix_get(r, 0, 0));

  for (k = 0 ; k < n ; k++)
    {
      double rkk = gsl_matrix_get(r, k, k);

      if (fabs(rkk) <= tolr)
        {
          break;
        }

      gsl_matrix_set(r, k, k, 1.0/rkk);

      for (j = 0; j < k ; j++)
        {
          double t = gsl_matrix_get(r, j, k) / rkk;
          gsl_matrix_set (r, j, k, 0.0);

          for (i = 0; i <= j; i++)
            {
              double rik = gsl_matrix_get (r, i, k);
              double rij = gsl_matrix_get (r, i, j);

              gsl_matrix_set (r, i, k, rik - t * rij);
            }
        }
      kmax = k;
    }

  /* Form the full upper triangle of the inverse of R^T R in the full
     upper triangle of R */

  for (k = 0; k <= kmax ; k++)
    {
      for (j = 0; j < k; j++)
        {
          double rjk = gsl_matrix_get (r, j, k);

          for (i = 0; i <= j ; i++)
            {
              double rij = gsl_matrix_get (r, i, j);
              double rik = gsl_matrix_get (r, i, k);

              gsl_matrix_set (r, i, j, rij + rjk * rik);
            }
        }

      {
        double t = gsl_matrix_get (r, k, k);

        for (i = 0; i <= k; i++)
          {
            double rik = gsl_matrix_get (r, i, k);

            gsl_matrix_set (r, i, k, t * rik);
          };
      }
    }

  /* Form the full lower triangle of the covariance matrix in the
     strict lower triangle of R and in w */

  for (j = 0 ; j < n ; j++)
    {
      size_t pj = gsl_permutation_get (perm, j);

      for (i = 0; i <= j; i++)
        {
          size_t pi = gsl_permutation_get (perm, i);

          double rij;

          if (j > kmax)
            {
              gsl_matrix_set (r, i, j, 0.0);
              rij = 0.0 ;
            }
          else
            {
              rij = gsl_matrix_get (r, i, j);
            }

          if (pi > pj)
            {
              gsl_matrix_set (r, pi, pj, rij);
            }
          else if (pi < pj)
            {
              gsl_matrix_set (r, pj, pi, rij);
            }

        }

      {
        double rjj = gsl_matrix_get (r, j, j);
        gsl_matrix_set (covar, pj, pj, rjj);
      }
    }

  /* symmetrize the covariance matrix */

  for (j = 0 ; j < n ; j++)
    {
      for (i = 0; i < j ; i++)
        {
          double rji = gsl_matrix_get (r, j, i);

          gsl_matrix_set (covar, j, i, rji);
          gsl_matrix_set (covar, i, j, rji);
        }
    }

  gsl_matrix_free (r);
  gsl_permutation_free (perm);
  gsl_vector_free (tau);
  gsl_vector_free (norm);

  return GSL_SUCCESS;
}

max_stars_repo_name: ruslankuzmin/julia


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2016 Couchbase, Inc
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */
/** \file
 * This file is internal to the inner workings of
 * Phosphor and is not intended for public consumption.
 */

#pragma once

#include <array>
#include <iterator>
#include <memory>
#include <thread>
#include <type_traits>
#include <vector>
#include <functional>

#include <gsl_p/iterator.h>

#include "phosphor/platform/core.h"
#include "trace_event.h"

namespace phosphor {

#ifndef PHOSPHOR_CHUNK_PAGE_COUNT
#define PHOSPHOR_CHUNK_PAGE_COUNT 1
#endif

    /**
     * TraceChunk represents an array of TraceEvents
     *
     * The TraceChunk should be used from a single thread to
     * store various events.
     */
    class PHOSPHOR_API TraceChunk {
    public:
        static constexpr auto chunk_page_count = PHOSPHOR_CHUNK_PAGE_COUNT;
        static constexpr auto page_size = 4096;
        static constexpr auto array_offset = 64;
        static constexpr auto chunk_size =
            (((page_size * chunk_page_count) - array_offset) /
             sizeof(TraceEvent));
        using event_array = std::array<TraceEvent, chunk_size>;

        using const_iterator = event_array::const_iterator;

        /**
         * Constructor for a TraceChunk
         */
        TraceChunk() = default;

        /**
         * Reset the state of the TraceChunk
         *
         * This should be called before the TraceChunk is first used
         * as TraceChunk is a trivial type and requires initialisation.
         */
        void reset();

        /**
         * Used for adding TraceEvents to the chunk
         *
         * @return A reference to a TraceEvent to be replaced
         */
        TraceEvent& addEvent();

        /**
         * Used for reviewing TraceEvents in the chunk
         *
         * Valid indexes are from 0 to `count()`. There is no
         * bounds checking.
         *
         * @return A const reference to a TraceEvent in the chunk
         *         that can be used to review the event data
         */
        const TraceEvent& operator[](const int index) const;

        /**
         * Determine if the chunk is full
         *
         * @return true if the chunk is full (and should be replaced)
         *         or false otherwise.
         */
        bool isFull() const;

        /**
         * Determine up to which index of events is initialised
         *
         * @return The number of initialised events in the chunk
         */
        size_t count() const;

        /**
         * @return Const iterator to the start of the chunk
         */
        const_iterator begin() const;

        /**
         * @return Const iterator to the last initialised event in the chunk
         */
        const_iterator end() const;

    private:
        unsigned short next_free;
        event_array chunk;
    };

    // Forward decl
    class StatsCallback;

    /**
     * Abstract base-class for a buffer of TraceEvents
     *
     * The TraceBuffer loans out TraceChunks to individual
     * threads to reduce lock-contention on event logging.
     *
     * This class is *not* thread-safe and should only be directly
     * interacted either with the global lock owned by the TraceLog
     * or after tracing has been finished.
     *
     * A trace buffer can be iterated over using C++11 style range
     * iteration:
     *
     *     for(const auto& event : trace_buffer) {
     *         std::cout << event << std::endl;
     *     }
     *
     * Alternatively more complex iteration can be accomplished by
     * using the iterators returned by TraceBuffer::begin() and
     * TraceBuffer::end().
     *
     * Iteration should not be attempted while chunks are loaned out.
     */
    class TraceBuffer {
    public:
        /**
         * Virtual destructor to allow subclasses to be cleaned up
         * properly.
         */
        virtual ~TraceBuffer() = default;

        /**
         * Used for getting a TraceChunk to add events to
         *
         * @return A pointer to a TraceChunk to insert events into or
         *         nullptr if the buffer is full.
         */
        virtual TraceChunk* getChunk() = 0;

        /**
         * Used for returning a TraceChunk once full
         *
         * For some buffer implementations this *may* be a no-op
         * but for others which might reuse chunks this can be
         * used to only reuse chunks that have been finished with.
         *
         * @param chunk The chunk to be returned
         */
        virtual void returnChunk(TraceChunk& chunk) = 0;

        /**
         * Determine if there are no remaining chunks left to be
         * used
         *
         * @return true if there are no chunks left or false
         *         otherwise
         */
        virtual bool isFull() const = 0;

        /**
         * Callback for retrieving stats from the Buffer implementation
         *
         * Implementations MUST supply the following stats as minimum:
         *
         * - buffer_name <cstring_span>: Textual representation of buffer type
         * - buffer_is_full <bool>: True if the buffer is full
         * - buffer_chunk_count <size_t>: Chunks that are returned or loaned
         * - buffer_loaned_chunks <size_t>: Currently loaned chunks
         * - buffer_total_loaned <size_t>: Count of all chunks ever loaned
         * - buffer_size <size_t>: Max number of chunks that fit in the buffer
         * - buffer_generation <size_t>: Generation number of the buffer
         *
         * On a non-rotating buffer, if buffer_chunk_count is equal to
         * buffer-size then that must suggest the buffer is full and
         * there are no more chunks to be loaned. On a rotating buffer
         * it suggests that chunks are being reused.
         *
         * Buffer implementations may include other relevant stats but
         * end-users SHOULD NOT assume the existence of those stats.
         */
        virtual void getStats(StatsCallback& addStats) const = 0;

        /**
         * Used for accessing TraceChunks in the buffer
         *
         * Valid indexes are from 0 to `count()`. There is no
         * bounds checking.
         *
         * @return A const reference to a TraceEvent in the chunk
         *         that can be used to review the event data
         * @throw std::logic_error if chunks are currently loaned
         *        out to chunk tenants.
         */
        virtual const TraceChunk& operator[](const int index) const = 0;

        /**
         * Used for determining the number of chunks in the buffer
         *
         * @return Number of chunks currently in the buffer
         */
        virtual size_t chunk_count() const = 0;

        /**
         * @return The generation number of the TraceBuffer
         */
        virtual size_t getGeneration() const = 0;

        /**
         * Const bi-directional iterator over the TraceChunks in a TraceBuffer
         *
         * Usage:
         *
         *    chunk_iterator it = buffer.chunk_start();
         *     std::cout << *(it->start()) << std::endl;
         *
         *     for(const auto& chunk : buffer.chunks()) {
         *         for(const auto& event : chunk) {
         *             std::cout << event << std::endl;
         *         }
         *     }
         */
        class PHOSPHOR_API chunk_iterator
            : public std::iterator<std::bidirectional_iterator_tag,
                                   TraceChunk> {
            using const_reference = const TraceChunk&;
            using const_pointer = const TraceChunk*;
            using _self = chunk_iterator;

        public:
            chunk_iterator() = default;
            chunk_iterator(const TraceBuffer& buffer_);
            chunk_iterator(const TraceBuffer& buffer_, size_t index_);
            const_reference operator*() const;
            const_pointer operator->() const;
            chunk_iterator& operator++();
            bool operator==(const chunk_iterator& other) const;
            bool operator!=(const chunk_iterator& other) const;

        protected:
            const TraceBuffer& buffer;
            size_t index;
        };

        /**
         * @return A const iterator to the first chunk of the TraceBuffer
         */
        virtual chunk_iterator chunk_begin() const = 0;

        /**
         * @return A const iterator to after the last chunk of the TraceBuffer
         */
        virtual chunk_iterator chunk_end() const = 0;

        /**
         * Const iterator over a TraceBuffer, implements required methods of
         * a bi-directional iterator.
         *
         * Usage:
         *
         *     event_iterator it = trace_buffer.start();
         *     std::cout << *it << std::endl;
         *     std::cout << *(++it) << std::endl;
         *
         *     for(const auto& event : trace_buffer) {
         *         std::cout << event << std::endl;
         *     }
         *
         * For resource efficiency event_iterator does not have
         * post-increment
         */
        using event_iterator = gsl_p::multidimensional_iterator<chunk_iterator>;

        /**
         * @return A const iterator to the first event of the TraceBuffer
         */
        virtual event_iterator begin() const = 0;

        /**
         * @return A const iterator to after the last event of the TraceBuffer
         */
        virtual event_iterator end() const = 0;

        /**
         * Helper class that can be used to create for-range loops over
         * the chunks of the TraceBuffer
         *
         *     for(const auto& chunk : TraceBuffer::chunk_iterable(buffer) {
         *         // Do something with every chunk
         *     }
         */
        class PHOSPHOR_API chunk_iterable {
        public:
            /**
             * @param buffer_ The buffer to iterate over
             */
            chunk_iterable(const TraceBuffer& buffer_) : buffer(buffer_) {}

            /**
             * @return A const iterator to the first chunk of the buffer
             */
            chunk_iterator begin() {
                return buffer.chunk_begin();
            }

            /**
             * @return A const iterator to after the last chunk of the buffer
             */
            chunk_iterator end() {
                return buffer.chunk_end();
            }

        private:
            const TraceBuffer& buffer;
        };

        /**
         * Helper function for creating a chunk_iterable over the buffer
         *
         *     for(const auto& chunk : buffer.chunks() {
         *         // Do something with every chunk
         *     }
         *
         * @return An iterable over the class
         */
        chunk_iterable chunks() const {
            return chunk_iterable(*this);
        }
    };

    using buffer_ptr = std::unique_ptr<TraceBuffer>;

    /**
     * Interface for a TraceBuffer factory
     */
    using trace_buffer_factory =
        std::function<buffer_ptr(size_t generation, size_t buffer_size)>;

    PHOSPHOR_API
    buffer_ptr make_fixed_buffer(size_t generation, size_t buffer_size);

    PHOSPHOR_API
    buffer_ptr make_ring_buffer(size_t generation, size_t buffer_size);
}

ext: h

max_stars_repo_name: Chippiewill/Phosphor

max_issues_repo_name: Chippiewill/Phosphor

max_forks_repo_name: Chippiewill/Phosphor


text:
#include <stdio.h>
#include <math.h>
#include "mf.h"
#include <gsl/gsl_integration.h>

//****************************the functions are listed below**********
double func1 (double y, double T)  //The self-consistent equation f(y)=0
{
	double res = tanh(3*y/T)-y;
	return res;
}
double func2 (double y, double T) //The free energy per site
{
	double res = 0.5*3*y*y-T*log(2*cosh(3*y)/T);
	return res;
}
//*****************************************************************************
int main()
{
	FILE *fp;
	int anyi;
	if((fp=fopen(filename,"w"))==NULL)
	{
		printf("cannot open the file: %s\n",filename);
		printf("Please input an integer to end the program:");
		scanf("%d",&anyi);
		return 0;
	}

	FILE *ffp;
	if((ffp=fopen(filename2,"w"))==NULL)
	{
		printf("cannot open the file: %s\n",filename2);
		printf("Please input an integer to end the program:");
		scanf("%d",&anyi);
		return 0;
	}

	double tem0;          //judge indicator
	double T_step;             //the step of Temperature
	T_step = (T_high-T_low)/(N_T-1);    //note N_T is the tatal points of T
//	printf("The temperature step is %lf\n",T_step);

	double y_1;
	double y_2;

	double tem_1;
	double tem_2;

	double T = T_low;

	int i,j;    //loop indicator
	for(i=0;i<N_T;i++)
	{
		T = T_low + i*T_step;

		y_1=1e-1;
		y_2=1e7;
		tem_1 =  func1(y_1, T);
		tem_2 =  func1(y_2, T);
		//if(tem_1<0||tem_2>0)
		if(tem_1*tem_2>0)
		{
			printf("The two points in the boundary isn't good! Give up them now. Do it, baby!\n");
			return 0;
		}

		//******************************************

		for(j=10;j>1;j++)
		{
			tem_1 = func1(y_1, T);
			tem_2 = func1(y_2, T);
			tem0 = func1((y_1+y_2)/2.,T);
		//	printf("tem_1:%lf\ttem_2:%lf\n",tem_1,tem_2);
			if((tem_1<precision&&tem_1>=0)||(tem_1>-precision&&tem_1<=0))
			{
				fprintf(fp,"%lf\t%lf\n",T,y_1);
				fprintf(ffp,"%lf\t%lf\t%lf\n",T,func2(y_1,T),func2(-y_1,T));
				j=0;
			}
			if((tem_2<precision&&tem_2>=0)||(tem_2>-precision&&tem_2<=0))
			{
				fprintf(fp,"%lf\t%lf\n",T,y_2);
				fprintf(ffp,"%lf\t%lf\t%lf\n",T,func2(y_2,T),func2(-y_2,T));
				j=0;
			}
			if(tem_1*tem0>0)
			{
				y_1 = (y_1+y_2)/2;
			}
			if(tem_2*tem0>0)
			{
				y_2 = (y_1+y_2)/2;
			}
			printf("%lf\n",(y_1+y_2)/2); //the mediate point value

		}
	}
	return 1;
}

max_stars_repo_name: PostQuantum/Buckyball-Ising-Model

max_issues_repo_name: PostQuantum/Buckyball-Ising-Model

max_forks_repo_name: PostQuantum/Buckyball-Ising-Model


text:
/**
 *
 * @file core_ztsrfb.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Azzam Haidar
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <cblas.h>
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_ztsrfb applies a complex block reflector H or its transpose H' to a
 *  complex rectangular matrix formed by coupling two tiles A1 and A2,
 *  from either the left or the right.
 *
 *******************************************************************************
 *
 * @param[in] side
 *         @arg PlasmaLeft  : apply Q or Q**H from the Left;
 *         @arg PlasmaRight : apply Q or Q**H from the Right.
 *
 * @param[in] trans
 *         @arg PlasmaNoTrans   : No transpose, apply Q;
 *         @arg PlasmaConjTrans : ConjTranspose, apply Q**H.
 *
 * @param[in] direct
 *         Indicates how H is formed from a product of elementary
 *         reflectors
 *         @arg PlasmaForward  : H = H(1) H(2) . . . H(k) (Forward)
 *         @arg PlasmaBackward : H = H(k) . . . H(2) H(1) (Backward)
 *
 * @param[in] storev
 *         Indicates how the vectors which define the elementary
 *         reflectors are stored:
 *         @arg PlasmaColumnwise
 *         @arg PlasmaRowwise
 *
 * @param[in] M1
 *         The number of columns of the tile A1. M1 >= 0.
 *
 * @param[in] N1
 *         The number of rows of the tile A1. N1 >= 0.
 *
 * @param[in] M2
 *         The number of columns of the tile A2. M2 >= 0.
 *
 * @param[in] N2
 *         The number of rows of the tile A2. N2 >= 0.
 *
 * @param[in] K
 *          The order of the matrix T (= the number of elementary
 *          reflectors whose product defines the block reflector).
 *
 * @param[in,out] A1
 *         On entry, the M1-by-N1 tile A1.
 *         On exit, A1 is overwritten by the application of Q.
 *
 * @param[in] LDA1
 *         The leading dimension of the array A1. LDA1 >= max(1,N1).
 *
 * @param[in,out] A2
 *         On entry, the M2-by-N2 tile A2.
 *         On exit, A2 is overwritten by the application of Q.
 *
 * @param[in] LDA2
 *         The leading dimension of the tile A2. LDA2 >= max(1,N2).
 *
 * @param[in] V
 *         (LDV,K) if STOREV = 'C'
 *         (LDV,M2) if STOREV = 'R' and SIDE = 'L'
 *         (LDV,N2) if STOREV = 'R' and SIDE = 'R'
 *         The matrix V.
 *
 * @param[in] LDV
 *         The leading dimension of the array V.
 *         If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M2);
 *         if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N2);
 *         if STOREV = 'R', LDV >= K.
 *
 * @param[out] T
 *         The triangular K-by-K matrix T in the representation of the
 *         block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] LDT
 *         The leading dimension of the array T. LDT >= K.
 *
 * @param[in,out] WORK
 *
 * @param[in] LDWORK
 *         The dimension of the array WORK.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_ztsrfb = PCORE_ztsrfb
#define CORE_ztsrfb PCORE_ztsrfb
#endif
int
CORE_ztsrfb(PLASMA_enum side, PLASMA_enum trans, int direct, int storev,
            int M1, int N1, int M2, int N2, int K,
            PLASMA_Complex64_t *A1, int LDA1,
            PLASMA_Complex64_t *A2, int LDA2,
            const PLASMA_Complex64_t *V, int LDV,
            const PLASMA_Complex64_t *T, int LDT,
            PLASMA_Complex64_t *WORK, int LDWORK)
{
    static PLASMA_Complex64_t zone  =  1.0;
    static PLASMA_Complex64_t mzone = -1.0;

    int j;

    /* Check input arguments */
    if (M1 < 0) {
        coreblas_error(5, "Illegal value of M1");
        return -5;
    }
    if (N1 < 0) {
        coreblas_error(6, "Illegal value of N1");
        return -6;
    }
    if ( (M2 < 0) ||
         ( (M2 != M1) && (side == PlasmaRight) ) ){
        coreblas_error(7, "Illegal value of M2");
        return -7;
    }
    if ( (N2 < 0) ||
         ( (N2 != N1) && (side == PlasmaLeft) ) ){
        coreblas_error(8, "Illegal value of N2");
        return -8;
    }
    if (K < 0) {
        coreblas_error(9, "Illegal value of K");
        return -9;
    }

    /* Quick return */
    if ((M1 == 0) || (N1 == 0) || (M2 == 0) || (N2 == 0) || (K == 0))
        return PLASMA_SUCCESS;

    if (storev == PlasmaColumnwise) {
        if (direct == PlasmaForward) {
            if (side == PlasmaLeft) {
                /*
                 * B = A1 + V' * A2
                 */
                LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
                    lapack_const(PlasmaUpperLower),
                    K, N1,
                    A1, LDA1, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasConjTrans, CblasNoTrans,
                    K, N2, M2,
                    CBLAS_SADDR(zone), V, LDV,
                    A2, LDA2,
                    CBLAS_SADDR(zone), WORK, LDWORK);
                /*
                 * A2 = A2 - V*T*B -> B = T*B, A2 = A2 - V*B
                 */
                cblas_ztrmm(
                    CblasColMajor, CblasLeft, CblasUpper,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, K, N2,
                    CBLAS_SADDR(zone), T, LDT, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasNoTrans, CblasNoTrans,
                    M2, N2, K,
                    CBLAS_SADDR(mzone), V, LDV,
                    WORK, LDWORK,
                    CBLAS_SADDR(zone), A2, LDA2);
                /*
                 * A1 = A1 - B
                 */
                for(j = 0; j < N1; j++) {
                    cblas_zaxpy(
                        K, CBLAS_SADDR(mzone),
                        &WORK[LDWORK*j], 1,
                        &A1[LDA1*j], 1);
                }
            }
            /*
             * Columnwise / Forward / Right
             */
            else {
                /*
                 * B = A1 + A2 * V
                 */
                LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
                    lapack_const(PlasmaUpperLower),
                    M1, K,
                    A1, LDA1, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasNoTrans, CblasNoTrans,
                    M2, K, N2,
                    CBLAS_SADDR(zone), A2, LDA2,
                    V, LDV,
                    CBLAS_SADDR(zone), WORK, LDWORK);
                /*
                 * A2 = A2 - B*T*V' -> B = B*T, A2 = A2 - B*V'
                 */
                cblas_ztrmm(
                    CblasColMajor, CblasRight, CblasUpper,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, M1, K,
                    CBLAS_SADDR(zone), T, LDT, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasNoTrans, CblasConjTrans,
                    M2, N2, K,
                    CBLAS_SADDR(mzone), WORK, LDWORK,
                    V, LDV,
                    CBLAS_SADDR(zone), A2, LDA2);
                /*
                 * A1 = A1 - B
                 */
                for(j = 0; j < K; j++) {
                    cblas_zaxpy(
                        M1, CBLAS_SADDR(mzone),
                        &WORK[LDWORK*j], 1,
                        &A1[LDA1*j], 1);
                }
            }
        }
        else {
            coreblas_error(3, "Not implemented (ColMajor / Backward / Left or Right)");
            return PLASMA_ERR_NOT_SUPPORTED;
        }
    }
    else {
        if (direct == PlasmaForward) {
            /*
             * Rowwise / Forward / Left
             */
            if (side == PlasmaLeft) {
                /*
                 * B = A1 + V * A2
                 */
                LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
                    lapack_const(PlasmaUpperLower),
                    K, N1,
                    A1, LDA1, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasNoTrans, CblasNoTrans,
                    K, N2, M2,
                    CBLAS_SADDR(zone), V, LDV,
                    A2, LDA2,
                    CBLAS_SADDR(zone), WORK, LDWORK);
                /*
                 * A2 = A2 - V'*T*B -> B = T*B, A2 = A2 - V'*B
                 */
                cblas_ztrmm(
                    CblasColMajor, CblasLeft, CblasUpper,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, K, N2,
                    CBLAS_SADDR(zone), T, LDT, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasConjTrans, CblasNoTrans,
                    M2, N2, K,
                    CBLAS_SADDR(mzone), V, LDV,
                    WORK, LDWORK,
                    CBLAS_SADDR(zone), A2, LDA2);
                /*
                 * A1 = A1 - B
                 */
                for(j=0; j<N1; j++) {
                    cblas_zaxpy(
                        K, CBLAS_SADDR(mzone),
                        &WORK[LDWORK*j], 1,
                        &A1[LDA1*j], 1);
                }
            }
            /*
             * Rowwise / Forward / Right
             */
            else {
                /*
                 * B = A1 + A2 * V'
                 */
                LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
                    lapack_const(PlasmaUpperLower),
                    M1, K,
                    A1, LDA1, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasNoTrans, CblasConjTrans,
                    M2, K, N2,
                    CBLAS_SADDR(zone), A2, LDA2,
                    V, LDV,
                    CBLAS_SADDR(zone), WORK, LDWORK);
                /*
                 * A2 = A2 - B*T*V -> B = B*T, A2 = A2 - B*V'
                 */
                cblas_ztrmm(
                    CblasColMajor, CblasRight, CblasUpper,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, M1, K,
                    CBLAS_SADDR(zone), T, LDT, WORK, LDWORK);

                cblas_zgemm(
                    CblasColMajor, CblasNoTrans, CblasNoTrans,
                    M2, N2, K,
                    CBLAS_SADDR(mzone), WORK, LDWORK,
                    V, LDV,
                    CBLAS_SADDR(zone), A2, LDA2);
                /*
                 * A1 = A1 - B
                 */
                for(j = 0; j < K; j++) {
                    cblas_zaxpy(
                        M1, CBLAS_SADDR(mzone),
                        &WORK[LDWORK*j], 1,
                        &A1[LDA1*j], 1);
                }
            }
        }
        else {
            coreblas_error(3, "Not implemented (RowMajor / Backward / Left or Right)");
            return PLASMA_ERR_NOT_SUPPORTED;
        }
    }
    return PLASMA_SUCCESS;
}