elemop.h

//
// elemop.h
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <cljanss@limitpt.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit 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 Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
 
#ifndef _math_scmat_elemop_h
#define _math_scmat_elemop_h
 
#include <cmath>
#include <stdexcept>
#include <functional>
#include <util/state/state.h>
#include <util/state/statein.h>
#include <util/state/stateout.h>
#include <util/group/message.h>
#include <math/scmat/blkiter.h>
#include <math/scmat/predicate.h>
 
namespace sc {
 
class SCMatrixBlock;
class SCMatrixBlockIter;
class SCMatrixRectBlock;
class SCMatrixLTriBlock;
class SCMatrixDiagBlock;
class SCVectorSimpleBlock;
class SCMatrixRectSubBlock;
class SCMatrixLTriSubBlock;
class SCMatrixDiagSubBlock;
class SCVectorSimpleSubBlock;
 
class SCMatrix;
class SymmSCMatrix;
class DiagSCMatrix;
class SCVector;
 
struct SCMatrixIterationRanges {
  enum Values {AllElements = 0,
  Columns = 1,
  Rows = 2
  };
};
struct SCElement {
    typedef double value_type;
    SCElement() : irow(null.irow), icol(null.icol), value(null.value) {}
    SCElement(int ir, int ic, value_type v) : irow(ir), icol(ic), value(v) {}
    int irow;
    int icol;
    value_type value;
 
    static SCElement null;
 
    operator value_type() const {
      return value;
    }
};
 
template <class BinaryPredicate>
struct SCElementBinaryPredicateAdapter : public std::binary_function<SCElement, SCElement, bool> {
  bool operator()(const sc::SCElement& a, const sc::SCElement& b) {
    return BinaryPredicate()(a.value,b.value);
  }
};
 
class SCElementOp: public SavableState {
  public:
    SCElementOp();
    SCElementOp(StateIn&s): SavableState(s) {}
    virtual ~SCElementOp();
    virtual int has_collect();
    virtual void defer_collect(int);
    virtual void collect(const Ref<MessageGrp>&);
    virtual void collect(const Ref<SCElementOp>&);
    virtual int has_side_effects();
 
    virtual bool threadsafe() const;
 
    virtual bool cloneable() const;
 
    virtual Ref<SCElementOp> clone();
 
    virtual void process(SCMatrixBlockIter&) = 0;
 
    void process_base(SCMatrixBlock*block);
 
    virtual void process_spec_rect(SCMatrixRectBlock*);
    virtual void process_spec_ltri(SCMatrixLTriBlock*);
    virtual void process_spec_diag(SCMatrixDiagBlock*);
    virtual void process_spec_vsimp(SCVectorSimpleBlock*);
    virtual void process_spec_rectsub(SCMatrixRectSubBlock*);
    virtual void process_spec_ltrisub(SCMatrixLTriSubBlock*);
    virtual void process_spec_diagsub(SCMatrixDiagSubBlock*);
    virtual void process_spec_vsimpsub(SCVectorSimpleSubBlock*);
};
 
class SCElementOp2: public SavableState {
  public:
    SCElementOp2();
    SCElementOp2(StateIn&s): SavableState(s) {}
    virtual ~SCElementOp2();
    virtual int has_collect();
    virtual void defer_collect(int);
    virtual int has_side_effects();
    virtual int has_side_effects_in_arg();
    virtual void collect(const Ref<MessageGrp>&);
    virtual void process(SCMatrixBlockIter&,SCMatrixBlockIter&) = 0;
    void process_base(SCMatrixBlock*,SCMatrixBlock*);
    virtual void process_spec_rect(SCMatrixRectBlock*,SCMatrixRectBlock*);
    virtual void process_spec_ltri(SCMatrixLTriBlock*,SCMatrixLTriBlock*);
    virtual void process_spec_diag(SCMatrixDiagBlock*,SCMatrixDiagBlock*);
    virtual void process_spec_vsimp(SCVectorSimpleBlock*,SCVectorSimpleBlock*);
};
 
class SCElementOp3: public SavableState {
  public:
    SCElementOp3();
    SCElementOp3(StateIn&s): SavableState(s) {}
    virtual ~SCElementOp3();
    virtual int has_collect();
    virtual void defer_collect(int);
    virtual int has_side_effects();
    virtual int has_side_effects_in_arg1();
    virtual int has_side_effects_in_arg2();
    virtual void collect(const Ref<MessageGrp>&);
    virtual void process(SCMatrixBlockIter&,
                         SCMatrixBlockIter&,
                         SCMatrixBlockIter&) = 0;
    void process_base(SCMatrixBlock*,SCMatrixBlock*,SCMatrixBlock*);
    virtual void process_spec_rect(SCMatrixRectBlock*,
                                   SCMatrixRectBlock*,
                                   SCMatrixRectBlock*);
    virtual void process_spec_ltri(SCMatrixLTriBlock*,
                                   SCMatrixLTriBlock*,
                                   SCMatrixLTriBlock*);
    virtual void process_spec_diag(SCMatrixDiagBlock*,
                                   SCMatrixDiagBlock*,
                                   SCMatrixDiagBlock*);
    virtual void process_spec_vsimp(SCVectorSimpleBlock*,
                                    SCVectorSimpleBlock*,
                                    SCVectorSimpleBlock*);
};
 
class SCElementScalarProduct: public SCElementOp2 {
  private:
    int deferred_;
    double product;
  public:
    SCElementScalarProduct();
    SCElementScalarProduct(StateIn&);
    ~SCElementScalarProduct();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&,SCMatrixBlockIter&);
    int has_collect();
    void defer_collect(int);
    void collect(const Ref<MessageGrp>&);
    double result();
    void init() { product = 0.0; }
};
 
 
class SCElementDestructiveProduct: public SCElementOp2 {
  public:
    SCElementDestructiveProduct();
    SCElementDestructiveProduct(StateIn&);
    ~SCElementDestructiveProduct();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&,SCMatrixBlockIter&);
};
 
class SCElementDAXPY: public SCElementOp2 {
  public:
    SCElementDAXPY(double a);
    SCElementDAXPY(StateIn&);
    ~SCElementDAXPY();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&,SCMatrixBlockIter&);
  private:
    double a_;
};
 
class SCElementScale: public SCElementOp {
  private:
    double scale;
  public:
    SCElementScale(double a);
    SCElementScale(StateIn&);
    ~SCElementScale();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
};
 
class SCElementRandomize: public SCElementOp {
  private:
    double assign;
  public:
    SCElementRandomize();
    SCElementRandomize(StateIn&);
    ~SCElementRandomize();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
};
 
class SCElementAssign: public SCElementOp {
  private:
    double assign;
  public:
    SCElementAssign(double a);
    SCElementAssign(StateIn&);
    ~SCElementAssign();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
};
 
class SCElementSquareRoot: public SCElementOp {
  public:
    SCElementSquareRoot();
    SCElementSquareRoot(double a);
    SCElementSquareRoot(StateIn&);
    ~SCElementSquareRoot();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
};
 
class SCElementInvert: public SCElementOp {
  private:
    double threshold_;
    int nbelowthreshold_;
    int deferred_;
  public:
    SCElementInvert(double threshold = 0.0);
    SCElementInvert(StateIn&);
    ~SCElementInvert();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
    int has_collect();
    void defer_collect(int);
    void collect(const Ref<MessageGrp>&);
    void collect(const Ref<SCElementOp>&);
    int result() { return nbelowthreshold_; }
};
 
 
class SCElementScaleDiagonal: public SCElementOp {
  private:
    double scale_diagonal;
  public:
    SCElementScaleDiagonal(double a);
    SCElementScaleDiagonal(StateIn&);
    ~SCElementScaleDiagonal();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
};
 
class SCElementShiftDiagonal: public SCElementOp {
  private:
    double shift_diagonal;
  public:
    SCElementShiftDiagonal(double a);
    SCElementShiftDiagonal(StateIn&);
    ~SCElementShiftDiagonal();
    int has_side_effects();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
};
 
template <class BinaryPredicate = std::less<double>, SCMatrixIterationRanges::Values IterationRanges = SCMatrixIterationRanges::AllElements>
class SCElementFindExtremum: public SCElementOp {
  public:
    typedef SCElementFindExtremum<BinaryPredicate, IterationRanges> this_type;
 
    SCElementFindExtremum(double init_value = 0.0):deferred_(0), init_value_(init_value), r(1, SCElement(0,0,init_value)) {}
    SCElementFindExtremum(StateIn& s): SCElementOp(s)
    {
      {
        int nchar; s.get(nchar);
        r.resize(nchar / sizeof(SCElement));
        char* r_char;
        s.get(r_char);
        s.get(init_value_);
        memcpy(&r[0], r_char, nchar);
        delete[] r_char;
      }
      s.get(deferred_);
    }
 
    ~SCElementFindExtremum() {}
 
    void save_data_state(StateOut& s)
    {
      {
        int nchar = r.size() * sizeof(SCElement);
        s.put(nchar);
        s.put(init_value_);
        s.put((char*)(&r[0]), nchar);
      }
      s.put(deferred_);
    }
 
    void process(SCMatrixBlockIter& i) {
      for (i.reset(); i; ++i) {
        int range;
        switch(IterationRanges) {
          case SCMatrixIterationRanges::AllElements: range = 0; break;
          case SCMatrixIterationRanges::Columns: range = i.j(); break;
          case SCMatrixIterationRanges::Rows: range = i.i(); break;
        }
        if (range >= r.size()) {
          const size_t old_size = r.size();
          r.resize(range+1);
          for(size_t i=old_size; i<r.size(); ++i) {
            r[i] = SCElement(0,0,init_value_);
          }
        }
        double value = i.get();
        if (Op(r[range].value, value)) {
          r[range] = SCElement(i.i(), i.j(), i.get());
        }
      }
    }
 
    int has_collect() { return 1; }
    void defer_collect(int h) { deferred_=h; }
 
    void collect(const Ref<MessageGrp>& msg) {
      if (!deferred_) {
        GrpCompareReduce<SCElement, SCElementBinaryPredicateAdapter<BinaryPredicate> > reduce;
        msg->reduce<SCElement>(&r[0], r.size(), reduce);
      }
    }
    void collect(const Ref<SCElementOp>& op) {
      throw std::runtime_error(
          "SCElementMaxAbs::collect(const Ref<SCElementOp> &): not implemented");
    }
 
    const std::vector<SCElement>& result() { return r; }
 
  private:
    int deferred_;
    double init_value_;
    std::vector<SCElement> r;
    BinaryPredicate Op;
 
    static ClassDesc class_desc_;
};
 
template <class BinaryPredicate, SCMatrixIterationRanges::Values IterationRanges>
ClassDesc SCElementFindExtremum<BinaryPredicate,IterationRanges>::class_desc_(
  typeid(this_type),
  typeid(this_type).name(),
  1,"public SCElementOp",
  0, 0, sc::create<typename SCElementFindExtremum<BinaryPredicate,IterationRanges>::this_type>);
 
class SCElementMaxAbs: public SCElementOp {
  private:
    int deferred_;
    double r;
  public:
    SCElementMaxAbs();
    SCElementMaxAbs(StateIn&);
    ~SCElementMaxAbs();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
    int has_collect();
    void defer_collect(int);
    void collect(const Ref<MessageGrp>&);
    void collect(const Ref<SCElementOp>&);
    double result();
};
 
 
class SCElementMinAbs: public SCElementOp {
  private:
    int deferred_;
    double r;
  public:
    // rinit must be greater than the magnitude of the smallest element
    SCElementMinAbs(double rinit);
    SCElementMinAbs(StateIn&);
    ~SCElementMinAbs();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
    int has_collect();
    void defer_collect(int);
    void collect(const Ref<MessageGrp>&);
    void collect(const Ref<SCElementOp>&);
    double result();
};
 
class SCElementSum: public SCElementOp {
  private:
    int deferred_;
    double r;
  public:
    SCElementSum();
    SCElementSum(StateIn&);
    ~SCElementSum();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
    int has_collect();
    void defer_collect(int);
    void collect(const Ref<MessageGrp>&);
    void collect(const Ref<SCElementOp>&);
    double result();
    void init() { r = 0.0; }
};
 
class SCElementKNorm: public SCElementOp {
  private:
    int deferred_;
    double r_;  // result
    unsigned int k_;  // norm parameter
 
    static double _process(SCMatrixBlockIter& i, int k);
    static double _process1(SCMatrixBlockIter& i);
    static double _process2(SCMatrixBlockIter& i);
 
  public:
    SCElementKNorm(unsigned int k = 2);
    SCElementKNorm(StateIn&);
    ~SCElementKNorm();
    void save_data_state(StateOut&);
    void process(SCMatrixBlockIter&);
    int has_collect();
    void defer_collect(int);
    void collect(const Ref<MessageGrp>&);
    void collect(const Ref<SCElementOp>&);
    double result();
    void init() { r_ = 0.0; }
};
 
class SCElementDot: public SCElementOp {
  private:
    double** avects;
    double** bvects;
    int length;
  public:
    SCElementDot(StateIn&);
    void save_data_state(StateOut&);
    SCElementDot(double**a, double**b, int length);
    void process(SCMatrixBlockIter&);
    int has_side_effects();
};
 
class SCElementAccumulateSCMatrix: public SCElementOp {
  private:
    SCMatrix *m;
  public:
    SCElementAccumulateSCMatrix(SCMatrix *);
    int has_side_effects();
    void process(SCMatrixBlockIter&);
};
 
class SCElementAccumulateSymmSCMatrix: public SCElementOp {
  private:
    SymmSCMatrix *m;
  public:
    SCElementAccumulateSymmSCMatrix(SymmSCMatrix *);
    int has_side_effects();
    void process(SCMatrixBlockIter&);
};
 
class SCElementAccumulateDiagSCMatrix: public SCElementOp {
  private:
    DiagSCMatrix *m;
  public:
    SCElementAccumulateDiagSCMatrix(DiagSCMatrix *);
    int has_side_effects();
    void process(SCMatrixBlockIter&);
};
 
class SCElementAccumulateSCVector: public SCElementOp {
  private:
    SCVector *m;
  public:
    SCElementAccumulateSCVector(SCVector *);
    int has_side_effects();
    void process(SCMatrixBlockIter&);
};
 
}
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: