home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2400 / realmat.cc < prev   
Encoding:
C/C++ Source or Header  |  1990-12-28  |  9.7 KB  |  366 lines
  1. /* ----------------------------------------------------------------------------
  2.  
  3. realmat: realmat.cc
  4.  
  5. realmat is a C++ matrix class. The header file realmat.h describes its use.
  6.  
  7. Copyright (C) 1990.
  8.  
  9. A. Ronald Gallant
  10. P.O. Box 5513 
  11. Raleigh NC 27650-5513 
  12. USA   
  13.  
  14. Permission to use, copy, modify, and distribute this software and its 
  15. documentation for any purpose and without fee is hereby granted, provided 
  16. that the above copyright notice appear in all copies and that both that 
  17. copyright notice and this permission notice appear in supporting 
  18. documentation.  
  19.  
  20. This software is provided "as is" without any expressed or implied warranty.
  21.  
  22. ---------------------------------------------------------------------------- */
  23. #include "realmat.h"
  24.  
  25. realmat::realmat(INTEGER r, INTEGER c, REAL* a)
  26. {
  27.   if (r<=0) error("Error, realmat::realmat, Number of rows not positive");
  28.   if (c<=0) error("Error, realmat::realmat, Number of columns not positive");
  29.   rows=r; 
  30.   cols=c; 
  31.   len=r*c; 
  32.   x = a;
  33. }
  34.  
  35.  
  36. void realmat::realmat_constructor(INTEGER r, INTEGER c, REAL fill_value)
  37. {
  38.   if (r<=0) error("Error, realmat::realmat, Number of rows not positive");
  39.   if (c<=0) error("Error, realmat::realmat, Number of columns not positive");
  40.   rows=r; 
  41.   cols=c; 
  42.   len=r*c; 
  43.   x = new REAL[len];
  44.   if (x == 0) error("Error, realmat::realmat, Operator new failed");
  45.   REAL* top = &(x[len-1]);
  46.   REAL* t = x;
  47.   while (t <= top) *t++ = fill_value;
  48. }
  49.  
  50. realmat::realmat(realmat& a)
  51. {
  52.     rows=a.rows;
  53.     cols=a.cols;
  54.     if(rows<=0) 
  55.       error("Error, realmat::realmat, Number of rows not positive");
  56.     if(cols<=0) 
  57.       error("Error, realmat::realmat, Number of columns not positive");
  58.     len=a.rows*a.cols; 
  59.     x = new REAL[len];
  60.     if (x == 0) error("Error, realmat::realmat, Operator new failed");
  61.     REAL* top = &(x[len-1]);
  62.     REAL* t = x;
  63.     REAL* u = a.x;
  64.     while (t <= top) *t++ = *u++;
  65. }
  66.  
  67. void realmat::resize_constructor(INTEGER r, INTEGER c, REAL fill_value)
  68. {
  69.   if (r<=0) error("Error, realmat::resize, Number of rows not positive");
  70.   if (c<=0) error("Error, realmat::resize, Number of columns not positive");
  71.   delete x;
  72.   rows=r; 
  73.   cols=c; 
  74.   len=r*c; 
  75.   x = new REAL[len];
  76.   if (x == 0) error("Error, realmat::resize, Operator new failed");
  77.   REAL* top = &(x[len-1]);
  78.   REAL* t = x;
  79.   while (t <= top) *t++ = fill_value;
  80. }
  81.  
  82. realmat& realmat::operator=(realmat& a)
  83. {
  84.   if (this != &a) {
  85.     delete x;
  86.     rows=a.rows;
  87.     cols=a.cols;
  88.     len=a.len;
  89.     x = new REAL[len];
  90.     if (x == 0) error("Error, realmat::operator=, Operator new failed");
  91.     REAL* top = &(x[len-1]);
  92.     REAL* t = x;
  93.     REAL* u = a.x;
  94.     while (t <= top) *t++ = *u++;
  95.   }
  96.   return *this;
  97. }
  98.  
  99. REAL& realmat::check1(INTEGER i)
  100. {
  101.   if ((1<=i) && (i<=len))
  102.     return x[i-1];
  103.   else
  104.     error("Error, realmat::check1, Index out of range");
  105.     exit(1); //This keeps the compiler from complaining;
  106. }
  107.  
  108. REAL& realmat::check2(INTEGER i, INTEGER j)
  109. {
  110.   if ((1<=i) && (1<=j) && (i<=rows) && (j<=cols))
  111.     return x[i + rows*j - rows - 1];  //  return x[rows*(j-1)+i-1]
  112.   else
  113.     error("Error, realmat::check2, Index out of range");
  114.     exit(1); //This keeps the compiler from complaining;
  115. }
  116.  
  117. realmat T(realmat& a)
  118. {
  119.   INTEGER newrows = a.cols;
  120.   INTEGER newcols = a.rows;
  121.   INTEGER newlen = newrows*newcols;
  122.   REAL*   newx = new REAL[newlen];
  123.   if (newx == 0) a.error("Error, realmat, T, Operator new failed");
  124.   INTEGER i,j;
  125.   for (j = 0; j < a.cols; j++)
  126.   for (i = 0; i < a.rows; i++)
  127.     newx[j + newrows*i] = a.x[i + a.rows*j];
  128.   realmat d(newrows,newcols,newx);
  129.   return d;
  130. }
  131.  
  132. realmat operator+(realmat& a, realmat& b)
  133. {
  134.   if ((a.cols == b.cols) && (a.rows == b.rows)) {
  135.     REAL* newx = new REAL[a.len];
  136.     if (newx == 0) a.error("Error, realmat, operator+, Operator new failed");
  137.     REAL* top = &(newx[a.len-1]);
  138.     REAL* t = newx;
  139.     REAL* u = a.x;
  140.     REAL* v = b.x;
  141.     while (t <= top) *t++ = *u++ + *v++;
  142.     realmat d(a.rows,a.cols,newx);
  143.     return d;
  144.   }
  145.   else
  146.     a.error("Error, realmat, operator+, Matrices not conformable."); 
  147.     exit(1); //This keeps the compiler from complaining;
  148. }
  149.  
  150. realmat operator+(realmat& a)
  151. {
  152.   REAL* newx = new REAL[a.len];
  153.   if (newx == 0) a.error("Error, realmat, operator+, Operator new failed");
  154.   REAL* top = &(newx[a.len-1]);
  155.   REAL* t = newx;
  156.   REAL* u = a.x;
  157.   while (t <= top) *t++ = *u++;
  158.   realmat d(a.rows,a.cols,newx);
  159.   return d;
  160. }
  161.  
  162. realmat operator-(realmat& a, realmat& b)
  163. {
  164.   if ((a.cols == b.cols) && (a.rows == b.rows)) {
  165.     REAL* newx = new REAL[a.len];
  166.     if (newx == 0) a.error("Error, realmat, operator-, Operator new failed");
  167.     REAL* top = &(newx[a.len-1]);
  168.     REAL* t = newx;
  169.     REAL* u = a.x;
  170.     REAL* v = b.x;
  171.     while (t <= top) *t++ = *u++ - *v++;
  172.     realmat d(a.rows,a.cols,newx);
  173.     return d;
  174.   }
  175.   else
  176.     a.error("Error, realmat, operator-, Matrices not conformable."); 
  177.     exit(1); //This keeps the compiler from complaining;
  178. }
  179.  
  180. realmat operator-(realmat& a)
  181. {
  182.   REAL* newx = new REAL[a.len];
  183.   if (newx == 0) a.error("Error, realmat, operator-, Operator new failed");
  184.   REAL* top = &(newx[a.len-1]);
  185.   REAL* t = newx;
  186.   REAL* u = a.x;
  187.   while (t <= top) *t++ = - *u++;
  188.   realmat d(a.rows,a.cols,newx);
  189.   return d;
  190. }
  191.  
  192. realmat operator*(realmat& a, realmat& b)
  193. {
  194.   if (a.cols == b.rows) {
  195.     INTEGER newlen = a.rows*b.cols;
  196.     REAL*   newx = new REAL[newlen];
  197.     if (newx == 0) a.error("Error, realmat, operator*, Operator new failed");
  198.     REAL    zero = 0.0;
  199.     REAL*   top = &(newx[newlen-1]);
  200.     REAL*   t = newx;
  201.     while (t <= top) *t++ = zero;
  202.     INTEGER i,j,k;
  203.     for (j = 0; j < b.cols; j++)
  204.     for (k = 0; k < a.cols; k++)
  205.     for (i = 0; i < a.rows; i++)
  206.       newx[i+a.rows*j] += a.x[i+a.rows*k] * b.x[k+b.rows*j];
  207.     realmat d(a.rows,b.cols,newx);
  208.     return d;
  209.   }
  210.   else {
  211.     a.error("Error, realmat, operator*, Matrices not conformable."); 
  212.     exit(1); //this keeps compiler from complaining
  213.   }
  214. }
  215.  
  216. realmat operator*(REAL& a, realmat& b)
  217. {
  218.   REAL*   newx = new REAL[b.len];
  219.   if (newx == 0) b.error("Error, realmat, operator*, Operator new failed");
  220.   REAL* top = &(newx[b.len-1]);
  221.   REAL* t = newx;
  222.   REAL* u = b.x;
  223.   while (t <= top) *t++ = a * *u++;
  224.   realmat d(b.rows,b.cols,newx);
  225.   return d;
  226. }
  227.  
  228. realmat operator*(INTEGER& a, realmat& b)
  229. {
  230.   REAL* newx = new REAL[b.len];
  231.   if (newx == 0) b.error("Error, realmat, operator*, Operator new failed");
  232.   REAL  f = a;
  233.   REAL* top = &(newx[b.len-1]);
  234.   REAL* t = newx;
  235.   REAL* u = b.x;
  236.   while (t <= top) *t++ = f * *u++;
  237.   realmat d(b.rows,b.cols,newx);
  238.   return d;
  239. }
  240.  
  241. void default_realmat_error_handler(const char* msg)
  242. {
  243.   cerr << msg << "\n";
  244.   exit(1);
  245. }
  246.  
  247. ONE_ARG_ERROR_HANDLER_T realmat_error_handler = default_realmat_error_handler;
  248.  
  249. ONE_ARG_ERROR_HANDLER_T set_realmat_error_handler(ONE_ARG_ERROR_HANDLER_T f)
  250. {
  251.   ONE_ARG_ERROR_HANDLER_T old = realmat_error_handler;
  252.   realmat_error_handler = f;
  253.   return old;
  254. }
  255.  
  256. void realmat::error(const char* msg)
  257. {
  258.   (*realmat_error_handler)(msg);
  259. }
  260.  
  261. ostream& operator<<(ostream& stream, realmat& a)
  262. {
  263.   char    line[256], eol[2] = {'\n','\0'}, bl[3] = {' ','\n','\0'};
  264.   char    *next,*save,*fcode;
  265.   REAL    f,af;
  266.   INTEGER i,j,linesize,maxcol,pad,start,stop,r,c;
  267.  
  268.   r=rows(a);
  269.   c=cols(a);
  270.  
  271.   ((LINESIZE<72)||(LINESIZE>133)) ? (linesize=133) : (linesize=LINESIZE);
  272.   maxcol = (linesize - 8)/12;
  273.   if (c < maxcol) maxcol = c;
  274.   pad = (linesize - 8 - 12*maxcol)/2 + 1;
  275.   memset(line,' ',pad);
  276.   save = line + pad;
  277.  
  278.   start=1;
  279.   do {
  280.     stop = start - 1 + maxcol;
  281.     if (stop > c) stop = c;
  282.  
  283.     stream << bl;
  284.     stream << bl;
  285.  
  286.     next = save;
  287.     memset(next,' ',6);
  288.     next = next + 6;
  289.  
  290.     for (j = start; j <= stop; j++) {
  291.       if (j < 1000)      {sprintf(next,"      Col%3i",j); next = next+12;}
  292.       else if (j <32760) {sprintf(next,"      C%5i"  ,j); next = next+12;}
  293.       else               {sprintf(next,"      TooBig"  ); next = next+12;}
  294.     }                                   
  295.  
  296.     memcpy(next,eol,2);
  297.                           
  298.     stream << line;
  299.     stream << bl;
  300.  
  301.     for (i = 1; i <= r; i++) {
  302.       next = save;
  303.       if (i < 1000)      {sprintf(next,"Row%3i",i); next = next+6;}
  304.       else if (i <32760) {sprintf(next,"R%5i"  ,i); next = next+6;}
  305.       else               {sprintf(next,"TooBig"  ); next = next+6;}
  306.       for (j = start; j <= stop; j++) {
  307.         fcode = "%12.3e";
  308.         f = a.elem(i,j);
  309.         af = f; if (f < 0.E0)  af = -f ;
  310.         if (af < 1.E+8 ) fcode = "%12.0f";
  311.         if (af < 1.E+5 ) fcode = "%12.1f";
  312.         if (af < 1.E+4 ) fcode = "%12.2f";
  313.         if (af < 1.E+3 ) fcode = "%12.3f";
  314.         if (af < 1.E+2 ) fcode = "%12.4f";
  315.         if (af < 1.E+1 ) fcode = "%12.5f";
  316.         if (af < 1.E+0 ) fcode = "%12.6f";
  317.         if (af < 1.E-1 ) fcode = "%12.7f";
  318.         if (af < 1.E-2 ) fcode = "%12.8f";
  319.         if (af < 1.E-4 ) fcode = "%12.3e";
  320.         if (af < 1.E-30) fcode = "%12.1f";
  321.         sprintf(next,fcode,f); next = next+12;
  322.       }
  323.       memcpy(next,eol,2);
  324.       stream << line;
  325.     }
  326.     start = stop + 1;
  327.   } while (stop < c);
  328.  
  329.   return stream;
  330. }
  331.  
  332. realmat invpsd(realmat& a, REAL eps)
  333. {
  334.   if (a.cols != a.rows) 
  335.     a.error("Error, realmat, invpsd, Matrix not square.");
  336.   INTEGER i;
  337.   for (i=0; i<a.rows; i++)
  338.     if (a.x[i+a.rows*i] < 0) 
  339.       a.error("Error, realmat, invpsd, Matrix not positive semi-definite.");
  340.  
  341.   INTEGER newrows = a.cols;
  342.   INTEGER newcols = a.rows;
  343.   INTEGER newlen = newrows*newcols;
  344.   REAL*   newx = new REAL[newlen];
  345.   if (newx == 0) a.error("Error, realmat, invpsd, Operator new failed");
  346.  
  347.   REAL* top = &(newx[a.len-1]);
  348.   REAL* t = newx;
  349.   REAL* u = a.x;
  350.     while (t <= top) *t++ = *u++;
  351.  
  352.   REAL* s = new REAL[newrows];
  353.   if (s == 0) a.error("Error, realmat, invpsd, Operator new failed");
  354.  
  355.   dcond(newx,newrows,s,0);
  356.   INTEGER ier = dsweep(newx,newrows,eps);
  357.   dcond(newx,newrows,s,1);
  358.  
  359.   ier++;  //This stops a compiler warning.
  360.  
  361.   delete s;
  362.  
  363.   realmat d(newrows,newcols,newx);
  364.   return d;
  365. }
  366.