home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / yamp / virtram.cpp < prev    next >
Encoding:
Text File  |  1993-01-11  |  15.4 KB  |  654 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20.  
  21.  
  22. // #include "virt.h"
  23.  
  24. //////////////////////////////////////////// string functions
  25. strtype::strtype(strtype &str)     // copy constructor
  26. { int len = MAXCHARS;
  27.     len = strlen(str.name);
  28.     len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
  29.     strncpy(name, str.name, len);
  30.     name[len] = '\0';
  31. }
  32. strtype::strtype(char *str)
  33. { int len = MAXCHARS;
  34.     len = strlen(str);
  35.     len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
  36.     strncpy(name, str, len);
  37.     name[len] = '\0';
  38. }
  39. strtype::strtype(void)
  40. {
  41.     name[0] = '\0';
  42. }
  43.  
  44. strtype strtype::operator + (strtype str)
  45. {
  46.     
  47.     int len1, len2;
  48.     len1 = strlen(name);
  49.     len2 = strlen(str.name);
  50.     int len =(len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 3 : len2;
  51.     len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
  52.     strncat(name, str.name, len);
  53.     name[len + len1] = '\0';
  54.     return *this;
  55. }
  56. strtype strtype::operator + (const char *str)
  57. {
  58.     int len1 = strlen(name);
  59.     int len2 = strlen(str);
  60.     int len = (len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 2 : len2;
  61.     len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
  62.     strncat(name, str, len);
  63.     name[len + len1] = '\0';
  64.     return *this;
  65. }
  66. strtype strtype::operator = (strtype str)
  67. {
  68.     int len = (strlen(str.name) >= MAXCHARS) ? MAXCHARS - 3
  69.             : strlen(str.name);
  70.     strncpy(name, str.name, len);
  71.     name[len] = '\0';
  72.     return *this;
  73. }
  74. strtype strtype::operator = (const char *str)
  75. {
  76.     int len = (strlen(str) >= MAXCHARS) ? MAXCHARS - 2 : strlen(str);
  77.     strncpy(name, str, len);
  78.     name[len] = '\0';
  79.     return *this;
  80. }
  81.  
  82.  
  83. //////////////////////////////////////////////////matrix stack functions
  84. MStack::MStack(void)
  85. {
  86.     static int cnter = 0;
  87.     next = NULL;
  88.     stackloc = 0;
  89.     level = 0;
  90.     if (!cnter) {
  91.         cnter = 1;
  92.         level = 1;
  93.     }
  94. }
  95.  
  96. void MStack::Inclevel(void)
  97. {
  98.     Dispatch->level++;
  99. }
  100.  
  101. void MStack::Declevel(void)
  102. {
  103.  
  104.     (Dispatch-> level)--;
  105.     if (Dispatch-> level < 0) {
  106.         printf("ERROR: LEVEL has been decremented too often\n");
  107.         exit(1);
  108.     }
  109. }
  110. void VMatrix::NewReference(VMatrix &ROp)
  111. {
  112.     signature = SIGNATURE;
  113.     r = ROp.r;
  114.     c = ROp.c;
  115.     // release the current header and share it with ROp.hdr
  116.     PurgeVectors();
  117.     v = ROp.v;
  118.     v->nrefs += 1;
  119.     mcheck = v->mm;
  120. }
  121.  
  122. void MStack::Push(VMatrix &ROp)
  123. {
  124.     ROp.Garbage("Push");
  125.     MStack *temp = new MStack;
  126.     temp->NewReference(ROp);
  127.     temp->Nameit(ROp.Getname(ROp));
  128.     temp->next = Dispatch-> next;
  129.     Dispatch->next = temp;
  130.     temp->level = Dispatch-> level;
  131.     temp->stackloc =++ (Dispatch->stackloc);
  132. }
  133.  
  134. VMatrix& MStack::Pop(void)
  135. {
  136.     Garbage("Pop");
  137.     if (Dispatch-> next && Dispatch-> stackloc) {
  138.         MStack *t = Dispatch->next;
  139.         Dispatch->next = Dispatch-> next-> next;
  140.         delete t;
  141.         (Dispatch-> stackloc)--;
  142.     }
  143.     else { Nrerror("POP: Stack is empty, cant pop any more\n"); }
  144.     return *Dispatch;
  145. }
  146.  
  147. void MStack::Cleanstack(void)
  148. {
  149.     while (Dispatch-> next-> level >= Dispatch-> level
  150.         && Dispatch->next-> next) Dispatch-> Pop();
  151.  
  152. }
  153.  
  154. void MStack::PrintStack(void)
  155. {
  156.     MStack *temp = Dispatch;
  157.     int t = 1 + Dispatch->stackloc;
  158.  
  159.     while (t--) {
  160.         printf("\n Stack matrix :location level r c:%d %d %d %d: name: ",
  161.             t, temp->level, temp-> r, temp-> c);
  162.         temp->Showname();
  163.         temp = temp->next;
  164.     }
  165.     printf("\n\n");
  166. }
  167.  
  168. //////////////////////////////////////////////////////matrix class
  169.  
  170. VMatrix::VMatrix(void)
  171. {
  172.     r = c = 1;
  173.     Nameit("t");
  174.     curvecind = 0;
  175.     signature = SIGNATURE;
  176.     SetupVectors(r, c);
  177. }
  178.  
  179. VMatrix::VMatrix(const char *str, int rr, int cc)
  180. {
  181.     r = rr;
  182.     c = cc;
  183.     Nameit(str);
  184.     curvecind = 0;
  185.     signature = SIGNATURE;
  186.     SetupVectors(r, c);
  187. }
  188. VMatrix::VMatrix(VMatrix &ROp)     // copy constructor
  189. {
  190.     ROp.Garbage("Copy Constructor");
  191.     r = ROp.r;
  192.     c = ROp.c;
  193.     name = ROp.name;
  194.     curvecind = 0;
  195.     signature = SIGNATURE;
  196.     SetupVectors(r, c);
  197.     for (int i = 1; i <= r;++ i) {
  198.         for (int j = 1; j <= c;++ j) M(i, j) = ROp.m(i, j);
  199.     }
  200.     Dispatch->Cleanstack();
  201. }
  202.  
  203. VMatrix::~VMatrix(void)
  204. {
  205.     PurgeVectors();
  206.     signature = 0;
  207.     r = 0;
  208.     c = 0;
  209. }
  210.  
  211. //////////////////////////////////// matrix member functions
  212.  
  213. double VMatrix::Nrerror(const char *errormsg)
  214. {
  215.     double x;
  216.     printf("\nMATRIX ERROR: %s \nexiting to system\n", errormsg);
  217.     x = 2;
  218.     exit(1);
  219.     return x;
  220. }
  221.  
  222. void VMatrix::Garbage(const char *errormsg)
  223. {
  224.     #ifndef NO_CHECKING
  225.     int errorint = 0;
  226.     // if( Dispatch->signature != SIGNATURE) errorint++;
  227.     if (signature != SIGNATURE) errorint++;
  228.     if (v-> mm[- 1] != SIGNATURE) {
  229.         printf("v->mm[-1]= %f\n", v->mm[- 1]);
  230.         errorint++;
  231.     }
  232.     if (v-> mm != mcheck) errorint++;
  233.     if (!v->mm) errorint++;
  234.     if (errorint) {
  235.         Dispatch->PrintStack();
  236.         printf("\nFunction name: %s", errormsg);
  237.         Nrerror("Operating on Garbage");
  238.     }
  239.     #endif
  240.     return;
  241. }
  242.  
  243. void VMatrix::SetupVectors(int rr, int cc)
  244. {
  245.  
  246.     unsigned int numele = (rr * cc) + 1;
  247.     unsigned int siz = sizeof (double);
  248.  
  249.     if (numele > 8190)
  250.         Nrerror("allocation failure 1, too many doubles");
  251.  
  252.     v = new vdoub;
  253.     v->nrefs = 1;
  254.  
  255.     if (!(v->mm = new double[numele]))
  256.         Nrerror("allocation failure 2 in SetupVector()");
  257.  
  258.     v->mm[0] = SIGNATURE;
  259.     (v-> mm)++;
  260.     mcheck = v->mm;
  261.     r = rr;
  262.     c = cc;
  263. }
  264.  
  265. void VMatrix::PurgeVectors(void)
  266. {
  267.     Garbage("PurgeVectors");
  268.  
  269.     v->nrefs -= 1;
  270.     if (v-> nrefs > 0)
  271.         return;
  272.  
  273.     (v-> mm)--;
  274.     if (v-> mm[0] == SIGNATURE) {
  275.         delete v->mm;
  276.     } else cout << "DID NOT DELETE mm!!!!!!!!! " << endl;
  277.     delete v;
  278.  
  279. }
  280.  
  281. void VMatrix::DisplayMat(void)
  282. {
  283.     Garbage("DisplayMat");
  284.     int wid = Getwid();
  285.     int dec = Getdec();
  286.     Showname();
  287.     printf("\n\n");
  288.     for (int i = 1; i <= r;++ i) {
  289.         for (int j = 1; j <= c;++ j) {
  290.             printf("%*.*lf ", wid, dec, m(i, j));
  291.         }
  292.         printf("\n");
  293.     }
  294.     printf("\n");
  295. }
  296. void VMatrix::InfoMat(void)
  297. {
  298.     Garbage("InfoMat");
  299.     printf("\n");
  300.     Showname();
  301.     printf("\nr c: %d %d\n", r, c);
  302. }
  303.  
  304. void VMatrix::Replace(VMatrix &ROp)
  305. {
  306.     ROp.Garbage("Replace");
  307.     if (r != ROp.r || c != ROp.c) {
  308.         PurgeVectors();
  309.         SetupVectors(ROp.r, ROp.c);
  310.     }
  311.     name = ROp.name;
  312.     for (int i = 1; i <= r;++ i) {
  313.         for (int j = 1; j <= c;++ j) M(i, j) = ROp.m(i, j);
  314.     }
  315. }
  316.  
  317. void VMatrix::operator = (VMatrix &ROp)
  318. {
  319.     Garbage("operator =");
  320.     ROp.Garbage("operator =");
  321.     Replace(ROp);
  322.     Dispatch->Cleanstack();
  323. }
  324.  
  325. ////////////////////// end matrix member functions
  326.  
  327. /////////////////// freind functions of matrix class
  328.  
  329. ///------------- addition
  330.  
  331. VMatrix& operator + (VMatrix &LOp, VMatrix &ROp)
  332. {
  333.     LOp.Garbage("operator +");
  334.     ROp.Garbage("operator +");
  335.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  336.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  337.         for (int i = 1; i <= LOp.r;++ i) {
  338.             for (int j = 1; j <= ROp.c;++ j) {
  339.                 temp->M(i, j) = LOp.m(i, j) + ROp.m(i, j);
  340.             }
  341.         }
  342.         temp->name = temp-> name + LOp.name + "+" + ROp.name + ")";
  343.         Dispatch->Push(*temp);
  344.         delete temp;
  345.     }
  346.     else ROp.Nrerror("Mismatched Matrix sizes in addition function\n");
  347.     return Dispatch->ReturnMat();
  348. }
  349. VMatrix& operator +(double scalar, VMatrix &ROp)
  350. {
  351.     ROp.Garbage("operator s+M");
  352.     strtype string;
  353.     char buffer[MAXCHARS] = { '\0' };
  354.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  355.         for (int i = 1; i <= ROp.r;++ i) {
  356.             for (int j = 1; j <= ROp.c;++ j) {
  357.                 temp->M(i, j) = scalar + ROp.m(i, j);
  358.             }
  359.         }
  360.         gcvt(scalar, NDECS, buffer);
  361.         string = buffer;
  362.         temp->name = temp-> name + string + "+" + ROp.name + ")";
  363.         Dispatch->Push(*temp);
  364.         delete temp;
  365.         return Dispatch->ReturnMat();
  366. }
  367. VMatrix& operator +(VMatrix &ROp, double scalar)
  368. {
  369.     ROp.Garbage("operator M+s");
  370.     strtype string;
  371.     char buffer[MAXCHARS] = { '\0' };
  372.  
  373.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  374.         for (int i = 1; i <= ROp.r;++ i) {
  375.             for (int j = 1; j <= ROp.c;++ j) {
  376.                 temp->M(i, j) = scalar + ROp.m(i, j);
  377.             }
  378.         }
  379.         gcvt(scalar, NDECS, buffer);
  380.         string = buffer;
  381.         temp->name = temp-> name + ROp.name + "+" + string + ")";
  382.         Dispatch->Push(*temp);
  383.         delete temp;
  384.         return Dispatch->ReturnMat();
  385. }
  386.  
  387. ////-------------subtraction
  388.  
  389. VMatrix& operator -(VMatrix &LOp, VMatrix &ROp)
  390. {
  391.     LOp.Garbage("operator M-M");
  392.     ROp.Garbage("operator M-M");
  393.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  394.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  395.         for (int i = 1; i <= LOp.r;++ i) {
  396.             for (int j = 1; j <= LOp.c;++ j) {
  397.                 temp->M(i, j) = LOp.m(i, j) - ROp.m(i, j);
  398.             }
  399.         }
  400.         temp->name = temp-> name + LOp.name + "-" + ROp.name + ")";
  401.         Dispatch->Push(*temp);
  402.         delete temp;
  403.     }
  404.     else Dispatch->Nrerror(
  405.             "Mismatched VMatrix sizes in subtraction function\n");
  406.     return Dispatch->ReturnMat();
  407. }
  408. VMatrix& operator -(double scalar, VMatrix &ROp)
  409. {
  410.     ROp.Garbage("operator s-M");
  411.     strtype string;
  412.     char buffer[MAXCHARS] = { '\0' };
  413.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  414.         for (int i = 1; i <= ROp.r;++ i) {
  415.             for (int j = 1; j <= ROp.c;++ j) {
  416.                 temp->M(i, j) = scalar - ROp.m(i, j);
  417.             }
  418.         }
  419.         gcvt(scalar, NDECS, buffer);
  420.         string = buffer;
  421.         temp->name = temp-> name + string + "-" + ROp.name + ")";
  422.         Dispatch->Push(*temp);
  423.         delete temp;
  424.         return Dispatch->ReturnMat();
  425. }
  426. VMatrix& operator -(VMatrix &ROp, double scalar)
  427. {
  428.     ROp.Garbage("operator M-s");
  429.     strtype string;
  430.     char buffer[MAXCHARS] = { '\0' };
  431.  
  432.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  433.         for (int i = 1; i <= ROp.r;++ i) {
  434.             for (int j = 1; j <= ROp.c;++ j) {
  435.                 temp->M(i, j) = ROp.m(i, j) - scalar;
  436.             }
  437.         }
  438.         gcvt(scalar, NDECS, buffer);
  439.         string = buffer;
  440.         temp->name = temp-> name + ROp.name + "-" + string + ")";
  441.         Dispatch->Push(*temp);
  442.         delete temp;
  443.         return Dispatch->ReturnMat();
  444. }
  445. //------------- unary minus
  446. VMatrix& operator -(VMatrix &ROp)
  447. {
  448.     ROp.Garbage("operator -");
  449.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  450.     for (int i = 1; i <= ROp.r;++ i) {
  451.         for (int j = 1; j <= ROp.c;++ j) {
  452.             temp->M(i, j) = -ROp.m(i, j);
  453.         }
  454.     }
  455.     temp->name = temp-> name + "-" + ROp.name + ")";
  456.     Dispatch->Push(*temp);
  457.     delete temp;
  458.     return Dispatch->ReturnMat();
  459. }
  460.  
  461. //-------------- multiplication
  462.  
  463. VMatrix& operator* (VMatrix &LOp, VMatrix &ROp)
  464. {
  465.     long double sum = 0;
  466.     LOp.Garbage("operator M*M");
  467.     ROp.Garbage("operator M*M");
  468.     if (LOp.c == ROp.r) {
  469.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  470.         for (int i = 1; i <= LOp.r;++ i) {
  471.             for (int j = 1; j <= ROp.c;++ j) {
  472.                 sum = 0;
  473.                 for (int k = 1; k <= LOp.c; k++)
  474.                     sum +=(long double)
  475.                         ((long double) LOp.m(i, k)) *((long double) ROp.m(k,j));
  476.                 temp->M(i, j) = (double) sum;
  477.             }
  478.         }
  479.         temp->name = temp-> name + LOp.name + "*" + ROp.name + ")";
  480.         Dispatch->Push(*temp);
  481.         delete temp;
  482.     }
  483.     else Dispatch->Nrerror(
  484.             "Mismatched VMatrix sizes in multiplication function\n");
  485.  
  486.     return Dispatch->ReturnMat();
  487. }
  488. VMatrix& operator* (double scalar, VMatrix &ROp)
  489. {
  490.     ROp.Garbage("operator s*M");
  491.     strtype string;
  492.     char buffer[MAXCHARS] = { '\0' };
  493.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  494.         for (int i = 1; i <= ROp.r;++ i) {
  495.             for (int j = 1; j <= ROp.c;++ j) {
  496.                 temp->M(i, j) = scalar * ROp.m(i, j);
  497.             }
  498.         }
  499.         gcvt(scalar, NDECS, buffer);
  500.         string = buffer;
  501.         temp->name = temp-> name + string + "*" + ROp.name + ")";
  502.         Dispatch->Push(*temp);
  503.         delete temp;
  504.         return Dispatch->ReturnMat();
  505. }
  506. VMatrix& operator* (VMatrix &ROp, double scalar)
  507. {
  508.     ROp.Garbage("operator M*s");
  509.     strtype string;
  510.     char buffer[MAXCHARS] = { '\0' };
  511.  
  512.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  513.         for (int i = 1; i <= ROp.r;++ i) {
  514.             for (int j = 1; j <= ROp.c;++ j) {
  515.                 temp->M(i, j) = scalar * ROp.m(i, j);
  516.             }
  517.         }
  518.         gcvt(scalar, NDECS, buffer);
  519.         string = buffer;
  520.         temp->name = temp-> name + ROp.name + "*" + string + ")";
  521.         Dispatch->Push(*temp);
  522.         delete temp;
  523.         return Dispatch->ReturnMat();
  524. }
  525.  
  526. //-------- elementwise multiplication
  527.  
  528. VMatrix& operator %(VMatrix &LOp, VMatrix &ROp)
  529. {
  530.     LOp.Garbage("operator M%M");
  531.     ROp.Garbage("operator M%M");
  532.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  533.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  534.         for (int i = 1; i <= LOp.r;++ i) {
  535.             for (int j = 1; j <= ROp.c;++ j)
  536.                 temp->M(i, j) = LOp.m(i, j) *ROp.m(i, j);
  537.         }
  538.         temp->name = temp-> name + LOp.name + "%" + ROp.name + ")";
  539.         Dispatch->Push(*temp);
  540.         delete temp;
  541.     }
  542.     else Dispatch->Nrerror(
  543.             "Mismatched Matrix sizes in elementwise multiplication\n");
  544.     return Dispatch->ReturnMat();
  545. }
  546.  
  547. //------------- division
  548.  
  549. VMatrix& operator /(VMatrix &LOp, VMatrix &ROp)
  550. {
  551.     LOp.Garbage("operator M/M");
  552.     ROp.Garbage("operator M/M");
  553.  
  554.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  555.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  556.         for (int i = 1; i <= LOp.r;++ i) {
  557.             for (int j = 1; j <= ROp.c;++ j) {
  558.                 double d = ROp.m(i, j);
  559.                 double L = LOp.m(i, j);
  560.                 temp->M(i, j) = (fabs(d) > ZERO || fabs((d - L) / d) <     1.0E-5 )
  561.                               ? L / d
  562.                               :    ROp.Nrerror(" zero divisor in elementwise division");
  563.             }
  564.         }
  565.         temp->name = temp-> name + LOp.name + "/" + ROp.name + ")";
  566.         Dispatch->Push(*temp);
  567.         delete temp;
  568.     }
  569.     else Dispatch->Nrerror(
  570.             "Mismatched Matrix sizes in elementwise division\n");
  571.     return Dispatch->ReturnMat();
  572. }
  573.  
  574. VMatrix& operator /(VMatrix &ROp, double scalar)
  575. {
  576.     ROp.Garbage("operator M/s");
  577.     strtype string;
  578.     char buffer[MAXCHARS] = { '\0' };
  579.  
  580.         if (fabs(scalar) < ZERO)
  581.             ROp.Nrerror(" zero divisor in elementwise division");
  582.  
  583.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  584.         for (int i = 1; i <= ROp.r;++ i) {
  585.             for (int j = 1; j <= ROp.c;++ j) {
  586.                 temp->M(i, j) = ROp.m(i, j) / scalar;
  587.             }
  588.         }
  589.         gcvt(scalar, NDECS, buffer);
  590.         string = buffer;
  591.         temp->name = temp-> name + ROp.name + "/" + string + ")";
  592.         Dispatch->Push(*temp);
  593.         delete temp;
  594.         return Dispatch->ReturnMat();
  595. }
  596. VMatrix& operator /(double scalar, VMatrix &ROp)
  597. {
  598.     ROp.Garbage("operator s/M");
  599.     strtype string;
  600.     char buffer[MAXCHARS] = { '\0' };
  601.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  602.         for (int i = 1; i <= ROp.r;++ i) {
  603.             for (int j = 1; j <= ROp.c;++ j) {
  604.                 temp->M(i, j) = (ROp.m(i, j) != 0.0) ? scalar / ROp.m(i, j)
  605.                 : ROp.Nrerror("zero divisor in scalar divison");
  606.             }
  607.         }
  608.         gcvt(scalar, NDECS, buffer);
  609.         string = buffer;
  610.         temp->name = temp-> name + string + "/" + ROp.name + ")";
  611.         Dispatch->Push(*temp);
  612.         delete temp;
  613.         return Dispatch->ReturnMat();
  614. }
  615.  
  616. //////////////////////end of friend functions
  617.  
  618. ////////////////////// matrix functions
  619.  
  620. ////////////// utility functions
  621.  
  622. void VMatrix::Writeb(char *fid, VMatrix &mat)
  623.     /* write a matrix in an binary file */
  624.     {
  625.         int i;
  626.         FILE *file;
  627.         double x;
  628.         char name[MAXCHARS] = { '\0' };
  629.  
  630.         file = fopen(fid, "wb");
  631.         if (!file) Dispatch->Nrerror("WRITEB: unable to open file");
  632.  
  633.         mat.Garbage("Writeb");
  634.  
  635.         strncpy(name, mat.StringAddress(), 79);
  636.         i = strlen(name);
  637.         fwrite(&i, sizeof (int), 1, file);
  638.         fputs(name, file);
  639.         i = mat.r;
  640.         fwrite(&i, sizeof (int), 1, file);
  641.         i = mat.c;
  642.         fwrite(&i, sizeof (int), 1, file);
  643.  
  644.         for (i = 1; i <= mat.r; i++) {
  645.             for (int j = 1; j <= mat.c; j++){
  646.                 x = mat.m(i, j);
  647.                     fwrite(&x, sizeof (double), 1, file);
  648.                         }
  649.         }
  650.  
  651.         fclose(file);
  652.         mat.M(1, 1);     // reset matrix to base pointer
  653. }     /* writeb */
  654.