home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / nihcl-3.0 / ex / ArrayPartial.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-15  |  4.2 KB  |  177 lines
  1. // ArrayPartial.c  -- array of Partial automatic derivatives
  2.  
  3. // $Header: /afs/alw.nih.gov/unix/sun4_40c/usr/local/src/nihcl-3.0/share/ex/RCS/ArrayPartial.c,v 3.0 90/05/15 22:43:22 kgorlen Rel $
  4.  
  5. /*
  6. Author:
  7.     S. M. Orlow
  8.     Systex Inc.
  9.     Beltsville, MD 20705
  10.     301-4734-0111
  11.     sandy@alw.nih.gov
  12.  
  13. Function:
  14.  
  15. Class ArrayPartial is a container class used with class Partial. 
  16. The constructor, ArrayPartial::ArrayPartial(int size, double*) is 
  17. to be used only to construct the list of independent variables
  18. because it assumes that the number of variables is the same as the
  19. size of the list. 
  20. Typical usage of class ArrayPartial is:
  21.  
  22.     // NVAR = number of variables
  23.     #define NVAR 3
  24.     // NEQN = number of equations
  25.     #define NEQN 3
  26.     // x0 = vector of initial values
  27.     double x0[] = { 1.0, 1.0, 1.0 };
  28.         
  29.     // X[i] is the i-th variable
  30.     // (X[1][0],X[2][0],X[3][0]) is the solution vector
  31.     // initial value is (1,1,1)
  32.     ArrayPartial X(NVAR, x0);
  33.  
  34.     //f[i] is the i-th dependent variable
  35.     ArrayPartial f(NEQN);
  36. */
  37.  
  38. #include <libc.h>
  39.  
  40. #include <iostream.h>
  41.  
  42. #include "ArrayPartial.h"
  43. #include "Matrix.h"
  44.  
  45. // this constructor is only used to construct 
  46. // set of independent functions of one variable
  47. // assumes that siz is also the number of partials
  48. ArrayPartial::ArrayPartial(int siz, double* pd)
  49. {
  50.     sz = siz;
  51.     p = new Partial[sz];
  52.     for(int i=0; i<sz; i++) {
  53.       p[i].order(sz);
  54.       p[i][0] = (pd)? *pd:0;
  55.       p[i][i+1] = 1;  // i-th partial
  56.       }    
  57. }
  58. ArrayPartial::ArrayPartial(int siz, Partial* pp)
  59. {
  60.     sz = siz;
  61.     p = new Partial[sz];
  62.     if ( pp ) {
  63.       for(int i=0; i<sz; i++)
  64.         p[i] = pp[i];
  65.       }
  66.      else
  67.       for(int i=0; i<sz; i++)
  68.         p[i][i+1] = 1;
  69. }
  70. ArrayPartial::ArrayPartial(const ArrayPartial& ap)
  71. {
  72.     sz=ap.sz;
  73.     p = new Partial[sz];
  74.     for(int i=0; i<sz; i++)
  75.       p[i] = ap.at(i);
  76. }
  77. ArrayPartial::~ArrayPartial()
  78. {
  79.     delete p;
  80. }
  81.  
  82. void ArrayPartial::operator=(const ArrayPartial& ap)
  83. {
  84.     if ( sz !=ap.sz ) {
  85.       cerr << "operator=: ArrayPartial of size " << sz
  86.            << " expected, " << ap.sz << " found" << endl;
  87.       abort();
  88.       }
  89.     delete p;
  90.     p = ap.p;
  91. }
  92.  
  93. Partial& ArrayPartial::operator[](int Kth) const
  94. {
  95.     if ( Kth<1||Kth>sz ) {
  96.        cerr << "operator[]: index out of range" << endl;
  97.        abort();
  98.        }
  99.     return p[Kth-1];
  100. }
  101.  
  102. void ArrayPartial::operator+=(const Matrix& argm)
  103. {
  104.     Matrix m = argm;
  105.     if ( argm.nRow()==1 ) m = argm.t();
  106.     if ( m.nCol()!=1 ) {
  107.       cerr << "operator+=: matrix must be row or column vector"
  108.         << endl;
  109.       abort();
  110.       }
  111.     if ( m.nRow()!=sz ) {
  112.       cerr << "operator+=: vector length " << sz
  113.            << " expected, " << m.nRow() << " found" << endl;    
  114.       abort();
  115.       }
  116.     for(int i=0; i<sz; i++)
  117.       p[i][0] += m.at(i,0);
  118. }
  119. Matrix ArrayPartial::value() const
  120. {
  121.     Matrix val(sz,1);
  122.     for(int i=0; i<sz; i++)
  123.          val.at(i,0) = p[i][0];
  124.     return val;
  125. }
  126. Matrix ArrayPartial::jacobian() const
  127. {
  128.     Matrix val(sz,sz);
  129.     for(int i=0; i<sz; i++)
  130.       for(int j=0; j<sz; j++)
  131.          val.at(i,j) = p[i][j+1];
  132.     return val;
  133. }
  134. ArrayPartial ArrayPartial::operator+(const ArrayPartial& ap)
  135. {
  136.     if ( size() != ap.size() ) {
  137.       cerr << "operator+: ArrayPartial of size " << size()
  138.            << " expected, " << ap.size() << " found" << endl;
  139.       abort();
  140.       }
  141.     ArrayPartial val(sz);
  142.     for(int i=0; i<sz; i++)
  143.       val.at(i) = p[i] + ap.at(i);
  144.     return val;
  145. }
  146. ArrayPartial operator*(double d,const ArrayPartial& ap)
  147. {
  148.     ArrayPartial val(ap.size());
  149.     for(int i=0;i<ap.size(); i++)
  150.        val.at(i) = d*ap.at(i);
  151.     return val;
  152. }
  153. ArrayPartial operator*(const Matrix& m, const ArrayPartial& ap)
  154. {
  155.     if ( ap.size()!=m.nCol() ) {
  156.       cerr << "operator* mismatch: "
  157.            << m.nRow() << "x" << m.nCol() << " matrix and " 
  158.            << ap.size() << " ArrayPartial" << endl;
  159.       abort();
  160.       }
  161.     ArrayPartial val(ap.size());
  162.     for(int i=0; i< ap.size(); i++) {
  163.       Partial tmp;
  164.       for(int j=0; j<m.nCol(); j++ )
  165.         tmp = tmp + m.at(i,j)*ap.at(j);
  166.       val.at(i) = tmp;
  167.       }
  168.     return val;
  169. }
  170. void ArrayPartial::printOn(ostream& strm) const
  171. {
  172.     strm << "ArrayPartial[" << endl;
  173.     for(int i=0; i<sz; i++)
  174.       strm << p[i] << endl;
  175.     strm << "]" << endl;
  176. }
  177.