home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / math / euler / source / helpf.c < prev    next >
C/C++ Source or Header  |  1993-01-25  |  20KB  |  788 lines
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <math.h>
  4. #include <limits.h>
  5.  
  6. #include "header.h"
  7. #include "funcs.h"
  8. #include "helpf.h"
  9. #include "matheh.h"
  10.  
  11. #define wrong_arg() { error=26; output("Wrong argument\n"); return; }
  12.  
  13. void minmax (double *x, LONG n, double *min, double *max, 
  14.     int *imin, int *imax)
  15. /***** minmax
  16.     compute the total minimum and maximum of n double numbers.
  17. *****/
  18. {    LONG i;
  19.     if (n==0)
  20.     {    *min=0; *max=0; *imin=0; *imax=0; return; }
  21.     *min=*x; *max=*x; *imin=0; *imax=0; x++;
  22.     for (i=1; i<n; i++)
  23.     {    if (*x<*min) { *min=*x; *imin=(int)i; }
  24.         else if (*x>*max) { *max=*x; *imax=(int)i; }
  25.         x++;
  26.     }
  27. }
  28.  
  29. void transpose (header *hd)
  30. /***** transpose 
  31.     transpose a matrix
  32. *****/
  33. {    header *hd1,*st=hd;
  34.     double *m,*m1,*mh;
  35.     int c,r,i,j;
  36.     hd=getvalue(hd); if (error) return;
  37.     if (hd->type==s_matrix)
  38.     {    getmatrix(hd,&r,&c,&m);
  39.         hd1=new_matrix(c,r,""); if (error) return;
  40.         m1=matrixof(hd1);
  41.         for (i=0; i<r; i++)
  42.         {    mh=m1+i;
  43.             for (j=0; j<c; j++)
  44.             {    *mh=*m++; mh+=r;
  45.             }
  46.         }
  47.     }
  48.     else if (hd->type==s_cmatrix)
  49.     {    getmatrix(hd,&r,&c,&m);
  50.         hd1=new_cmatrix(c,r,""); if (error) return;
  51.         m1=matrixof(hd1);
  52.         for (i=0; i<r; i++)
  53.         {   mh=m1+(LONG)2*i;
  54.             for (j=0; j<c; j++)
  55.             {    *mh=*m++; *(mh+1)=*m++;
  56.                 mh+=(LONG)2*r;
  57.             }
  58.         }
  59.     }
  60.     else if (hd->type==s_real || hd->type==s_complex)
  61.     {    hd1=hd;
  62.     }
  63.     else
  64.     {    error=24; output("\' not defined for this value!\n");
  65.         return;
  66.     }
  67.     moveresult(st,hd1);
  68. }
  69.  
  70. void vectorize (header *init, header *step, header *end)
  71. {    double vinit,vstep,vend,*m;
  72.     int count;
  73.     header *result,*st=init;
  74.     init=getvalue(init); step=getvalue(step); end=getvalue(end);
  75.     if (error) return;
  76.     if (init->type!=s_real || step->type!=s_real || end->type!=s_real)
  77.     {    output("The : allows only real arguments!\n");
  78.         error=15; return;
  79.     }
  80.     vinit=*realof(init); vstep=*realof(step); vend=*realof(end);
  81.     if (vstep==0)
  82.     {    output("A step size of 0 is not allowed in : ");
  83.         error=401; return;
  84.     }
  85.     if (fabs(vend-vinit)/fabs(vstep)+1+epsilon>INT_MAX)
  86.     {    output1("A vector can only have %d elements",INT_MAX);
  87.         error=402; return;
  88.     }
  89.     count=(int)(floor(fabs(vend-vinit)/fabs(vstep)+1+epsilon));
  90.     if ((vend>vinit && vstep<0) || (vend<vinit && vstep>0))
  91.         count=0;
  92.     result=new_matrix(1,count,""); if (error) return;
  93.     m=matrixof(result);
  94.     while (count>=0)
  95.     {    *m++=vinit;
  96.         vinit+=vstep;
  97.         count--;
  98.     }
  99.     moveresult(st,result);
  100. }
  101.  
  102. void mfft (header *hd)
  103. {    header *st=hd,*result;
  104.     double *m,*mr;
  105.     int r,c;
  106.     hd=getvalue(hd);
  107.     if (error) return;
  108.     if (hd->type==s_real || hd->type==s_matrix)
  109.     {    make_complex(st); hd=st; }
  110.     getmatrix(hd,&r,&c,&m);
  111.     if (r!=1) wrong_arg();
  112.     result=new_cmatrix(1,c,"");
  113.     mr=matrixof(result);
  114.     memmove((char *)mr,(char *)m,(LONG)2*c*sizeof(double));
  115.     fft(mr,c,-1);
  116.     moveresult(st,result);
  117. }
  118.  
  119. void mifft (header *hd)
  120. {    header *st=hd,*result;
  121.     double *m,*mr;
  122.     int r,c;
  123.     hd=getvalue(hd);
  124.     if (error) return;
  125.     if (hd->type==s_real || hd->type==s_matrix)
  126.     {    make_complex(st); hd=st; }
  127.     getmatrix(hd,&r,&c,&m);
  128.     if (r!=1) wrong_arg();
  129.     result=new_cmatrix(1,c,"");
  130.     mr=matrixof(result);
  131.     memmove((char *)mr,(char *)m,(LONG)2*c*sizeof(double));
  132.     fft(mr,c,1);
  133.     moveresult(st,result);
  134. }
  135.  
  136. void mtridiag (header *hd)
  137. {    header *result,*st=hd,*result1;
  138.     double *m,*mr;
  139.     int r,c,*rows,i;
  140.     hd=getvalue(hd); if (error) return;
  141.     if (hd->type==s_matrix)
  142.     {    getmatrix(hd,&c,&r,&m);
  143.         if (c!=r || c==0) wrong_arg();
  144.         result=new_matrix(c,c,""); if (error) return;
  145.         result1=new_matrix(1,c,""); if (error) return;
  146.         mr=matrixof(result);
  147.         memmove(mr,m,(LONG)c*c*sizeof(double));
  148.         tridiag(mr,c,&rows);
  149.         mr=matrixof(result1);
  150.         for (i=0; i<c; i++) *mr++=rows[i]+1;
  151.     }
  152.     else if (hd->type==s_cmatrix)
  153.     {    getmatrix(hd,&c,&r,&m);
  154.         if (c!=r || c==0) wrong_arg();
  155.         result=new_cmatrix(c,c,""); if (error) return;
  156.         result1=new_matrix(1,c,""); if (error) return;
  157.         mr=matrixof(result);
  158.         memmove(mr,m,(LONG)c*c*(LONG)2*sizeof(double));
  159.         ctridiag(mr,c,&rows);
  160.         mr=matrixof(result1);
  161.         for (i=0; i<c; i++) *mr++=rows[i]+1;
  162.     }
  163.     else wrong_arg();
  164.     moveresult(st,result);
  165.     moveresult((header *)newram,result1);
  166. }
  167.  
  168. void mcharpoly (header *hd)
  169. {    header *result,*st=hd,*result1;
  170.     double *m,*mr;
  171.     int r,c;
  172.     hd=getvalue(hd); if (error) return;
  173.     if (hd->type==s_matrix)
  174.     {    getmatrix(hd,&c,&r,&m);
  175.         if (c!=r || c==0) wrong_arg();
  176.         result=new_matrix(c,c,""); if (error) return;
  177.         result1=new_matrix(1,c+1,""); if (error) return;
  178.         mr=matrixof(result);
  179.         memmove(mr,m,(LONG)c*c*sizeof(double));
  180.         charpoly(mr,c,matrixof(result1));
  181.     }
  182.     else if (hd->type==s_cmatrix)
  183.     {    getmatrix(hd,&c,&r,&m);
  184.         if (c!=r || c==0) wrong_arg();
  185.         result=new_cmatrix(c,c,""); if (error) return;
  186.         result1=new_cmatrix(1,c+1,""); if (error) return;
  187.         mr=matrixof(result);
  188.         memmove(mr,m,(LONG)c*c*(LONG)2*sizeof(double));
  189.         ccharpoly(mr,c,matrixof(result1));
  190.     }
  191.     else wrong_arg();
  192.     moveresult(st,result1);
  193. }
  194.  
  195. void mscompare (header *hd)
  196. {    header *st=hd,*hd1,*result;
  197.     hd=getvalue(hd);
  198.     hd1=getvalue(nextof(st));
  199.     if (error) return;
  200.     if (hd->type==s_string && hd1->type==s_string)
  201.     {    result=new_real(strcmp(stringof(hd),stringof(hd1)),"");
  202.         if (error) return;
  203.     }
  204.     else wrong_arg();
  205.     moveresult(st,result);
  206. }
  207.  
  208. void mfind (header *hd)
  209. {    header *st=hd,*hd1,*result;
  210.     double *m,*m1,*mr;
  211.     int i,j,k,c,r,c1,r1;
  212.     hd=getvalue(hd);
  213.     hd1=getvalue(nextof(st));
  214.     if (error) return;
  215.     if ((hd->type!=s_matrix && hd->type!=s_real) || 
  216.         (hd1->type!=s_matrix && hd1->type!=s_real)) wrong_arg();
  217.     getmatrix(hd,&c,&r,&m);
  218.     getmatrix(hd1,&c1,&r1,&m1);
  219.     if (c!=1 && r!=1) wrong_arg();
  220.     if (r!=1) c=r;
  221.     result=new_matrix(c1,r1,""); if (error) return;
  222.     mr=matrixof(result);
  223.     for (i=0; i<r1; i++)
  224.         for (j=0; j<c1; j++)
  225.         {    k=0;
  226.             while (k<c && m[k]<=*m1) k++;
  227.             if (k==c && *m1<=m[c-1]) k=c-1;
  228.             *mr++=k; m1++;
  229.         }
  230.     moveresult(st,result);
  231. }
  232.  
  233. void mdiag2 (header *hd)
  234. {    header *st=hd,*hd1,*result;
  235.     int c,r,i,n,l;
  236.     double *m,*mh,*mr;
  237.     hd1=next_param(hd);
  238.     hd=getvalue(hd); if (hd1) hd1=getvalue(hd1);
  239.     if (error) return;
  240.     if (hd1->type!=s_real) wrong_arg();
  241.     n=(int)*realof(hd1);
  242.     if (hd->type==s_matrix || hd->type==s_real)
  243.     {    getmatrix(hd,&r,&c,&m);
  244.         result=new_matrix(1,r,""); if (error) return;
  245.         mr=matrixof(result); l=0;
  246.         for (i=0; i<r; i++)
  247.         {    if (i+n>=c) break;
  248.             if (i+n>=0) { l++; *mr++=*mat(m,c,i,i+n); }
  249.         }
  250.         dimsof(result)->c=l;
  251.         result->size=matrixsize(1,c);
  252.     }
  253.     else if (hd->type==s_cmatrix || hd->type==s_complex)
  254.     {    getmatrix(hd,&r,&c,&m);
  255.         result=new_cmatrix(1,r,""); if (error) return;
  256.         mr=matrixof(result); l=0;
  257.         for (i=0; i<r; i++)
  258.         {    if (i+n>=c) break;
  259.             if (i+n>=0) 
  260.             {    l++;
  261.                 mh=cmat(m,c,i,i+n);
  262.                 *mr++=*mh++;
  263.                 *mr++=*mh;
  264.             }
  265.         }
  266.         dimsof(result)->c=l;
  267.         result->size=cmatrixsize(1,c);
  268.     }
  269.     else wrong_arg();
  270.     moveresult(st,result);
  271. }
  272.  
  273. void mband (header *hd)
  274. {    header *st=hd,*hd1,*hd2,*result;
  275.     int i,j,c,r,n1,n2;
  276.     double *m,*mr;
  277.     hd1=next_param(hd);
  278.     hd2=next_param(hd1);
  279.     hd=getvalue(hd); hd1=getvalue(hd1); hd2=getvalue(hd2);
  280.     if (error) return;
  281.     if (hd1->type!=s_real || hd2->type!=s_real) wrong_arg();
  282.     n1=(int)*realof(hd1); n2=(int)*realof(hd2);
  283.     if (hd->type==s_matrix || hd->type==s_real)
  284.     {    getmatrix(hd,&r,&c,&m);
  285.         result=new_matrix(r,c,""); if (error) return;
  286.         mr=matrixof(result);
  287.         for (i=0; i<r; i++)
  288.             for (j=0; j<c; j++)
  289.             {    if (j-i>=n1 && j-i<=n2) *mr++=*m++;
  290.                 else { *mr++=0.0; m++; }
  291.             }
  292.     }
  293.     else if (hd->type==s_cmatrix || hd->type==s_complex)
  294.     {    getmatrix(hd,&r,&c,&m);
  295.         result=new_cmatrix(r,c,""); if (error) return;
  296.         mr=matrixof(result);
  297.         for (i=0; i<r; i++)
  298.             for (j=0; j<c; j++)
  299.             {    if (j-i>=n1 && j-i<=n2) { *mr++=*m++; *mr++=*m++; }
  300.                 else { *mr++=0.0; *mr++=0.0; m+=2; }
  301.             }
  302.     }
  303.     else wrong_arg();
  304.     moveresult(st,result);
  305. }
  306.  
  307. void mdup (header *hd)
  308. {    header *result,*st=hd,*hd1;
  309.     double x,*m,*m1,*m2;
  310.     int c,i,n,j,r;
  311.     hd=getvalue(hd);
  312.     hd1=next_param(st);
  313.     if (!hd1) wrong_arg();
  314.     hd1=getvalue(hd1); if (error) return;
  315.     if (hd1->type!=s_real) wrong_arg();
  316.     x=*realof(hd1); n=(int)x;
  317.     if (n<1 || x>=INT_MAX) wrong_arg();
  318.     if (hd->type==s_matrix && dimsof(hd)->r==1)
  319.     {    c=dimsof(hd)->c;
  320.         result=new_matrix(n,c,"");
  321.         if (error) return;
  322.         m1=matrixof(hd); m2=matrixof(result);
  323.         for (i=0; i<n; i++)
  324.         {    m=mat(m2,c,i,0);
  325.             memmove((char *)m,(char *)m1,c*sizeof(double));
  326.         }
  327.     }
  328.     else if (hd->type==s_matrix && dimsof(hd)->c==1)
  329.     {    r=dimsof(hd)->r;
  330.         result=new_matrix(r,n,"");
  331.         if (error) return;
  332.         m1=matrixof(hd); m2=matrixof(result);
  333.         for (i=0; i<r; i++)
  334.         {    m=mat(m2,c,i,0);
  335.             for (j=0; j<n; j++)
  336.                 *mat(m2,n,i,j)=*mat(m1,1,i,0);
  337.         }
  338.     }
  339.     else if (hd->type==s_real)
  340.     {    result=new_matrix(n,1,""); if (error) return;
  341.         m1=matrixof(result);
  342.         for (i=0; i<n; i++) *m1++=*realof(hd);
  343.     }
  344.     else if (hd->type==s_cmatrix && dimsof(hd)->r==1)
  345.     {    c=dimsof(hd)->c;
  346.         result=new_cmatrix(n,c,"");
  347.         if (error) return;
  348.         m1=matrixof(hd); m2=matrixof(result);
  349.         for (i=0; i<n; i++)
  350.         {    m=cmat(m2,c,i,0);
  351.             memmove((char *)m,(char *)m1,(LONG)2*c*sizeof(double));
  352.         }
  353.     }
  354.     else if (hd->type==s_cmatrix && dimsof(hd)->c==1)
  355.     {    r=dimsof(hd)->r;
  356.         result=new_cmatrix(r,n,"");
  357.         if (error) return;
  358.         m1=matrixof(hd); m2=matrixof(result);
  359.         for (i=0; i<r; i++)
  360.         {    m=mat(m2,c,i,0);
  361.             for (j=0; j<n; j++)
  362.                 copy_complex(cmat(m2,n,i,j),cmat(m1,1,i,0));
  363.         }
  364.     }
  365.     else if (hd->type==s_complex)
  366.     {    result=new_cmatrix(n,1,""); if (error) return;
  367.         m1=matrixof(result);
  368.         for (i=0; i<n; i++) { *m1++=*realof(hd); *m1++=*imagof(hd); }
  369.     }
  370.     else wrong_arg();
  371.     moveresult(st,result);
  372. }
  373.  
  374. void make_complex (header *hd)
  375. /**** make_complex
  376.     make a function argument complex.
  377. ****/
  378. {    header *old=hd,*nextarg;
  379.     size_t size;
  380.     int r,c,i,j;
  381.     double *m,*m1;
  382.     hd=getvalue(hd);
  383.     if (hd->type==s_real)
  384.     {    size=sizeof(header)+2*sizeof(double);
  385.         nextarg=nextof(old);
  386.         if (newram+(size-old->size)>ramend)
  387.         {    output("Memory overflow!\n"); error=180; return; }
  388.         if (newram>(char *)nextarg)
  389.             memmove((char *)old+size,(char *)nextarg,
  390.                 newram-(char *)nextarg);
  391.         newram+=size-old->size;
  392.         *(old->name)=0; old->size=size;
  393.         old->type=s_complex;
  394.         *realof(old)=*realof(hd);
  395.         *imagof(old)=0.0;
  396.     }
  397.     else if (hd->type==s_matrix)
  398.     {    getmatrix(hd,&r,&c,&m);
  399.         size=cmatrixsize(r,c);
  400.         nextarg=nextof(old);
  401.         if (newram+(size-old->size)>ramend)
  402.         {    output("Memory overflow!\n"); error=180; return; }
  403.         if (newram>(char *)nextarg)
  404.             memmove((char *)old+size,(char *)nextarg,
  405.                 newram-(char *)nextarg);
  406.         newram+=size-old->size;
  407.         *(old->name)=0; old->size=size;
  408.         old->type=s_cmatrix;
  409.         dimsof(old)->r=r; dimsof(old)->c=c;
  410.         m1=matrixof(old);
  411.         for (i=r-1; i>=0; i--)
  412.             for (j=c-1; j>=0; j--)
  413.             {    *cmat(m1,c,i,j)=*mat(m,c,i,j);
  414.                 *(cmat(m1,c,i,j)+1)=0.0;
  415.             }
  416.     }
  417. }
  418.  
  419. void mvconcat (header *hd)
  420. {    header *st=hd,*hd1,*result;
  421.     double *m,*m1;
  422.     int r,c,r1,c1;
  423.     size_t size=0;
  424.     hd=getvalue(hd);
  425.     hd1=next_param(st);
  426.     if (!hd1) wrong_arg();
  427.     hd1=getvalue(hd1); if (error) return;
  428.     if (hd->type==s_real || hd->type==s_matrix)
  429.     {    if (hd1->type==s_real || hd1->type==s_matrix)
  430.         {    getmatrix(hd,&r,&c,&m);
  431.             getmatrix(hd1,&r1,&c1,&m1);
  432.             if (r!=0 && (c!=c1 || (LONG)r+r1>INT_MAX)) wrong_arg();
  433.             result=new_matrix(r+r1,c1,"");
  434.             if (r!=0)
  435.             {    size=(LONG)r*c*sizeof(double);
  436.                 memmove((char *)matrixof(result),(char *)m,
  437.                     size);
  438.             }
  439.             memmove((char *)matrixof(result)+size,
  440.                 (char *)m1,(LONG)r1*c1*sizeof(double));
  441.         }
  442.         else if (hd1->type==s_complex || hd1->type==s_cmatrix)
  443.         {    make_complex(st);
  444.             mvconcat(st);
  445.             return;
  446.         }
  447.         else wrong_arg();
  448.     }
  449.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  450.     {    if (hd1->type==s_complex || hd1->type==s_cmatrix)
  451.         {    getmatrix(hd,&r,&c,&m);
  452.             getmatrix(hd1,&r1,&c1,&m1);
  453.             if (r!=0 && (c!=c1 || (LONG)r+r1>INT_MAX)) wrong_arg();
  454.             result=new_cmatrix(r+r1,c1,"");
  455.             if (r!=0)
  456.             {    size=(LONG)r*(LONG)2*c*sizeof(double);
  457.                 memmove((char *)matrixof(result),(char *)m,
  458.                     size);
  459.             }
  460.             memmove((char *)matrixof(result)+size,
  461.                 (char *)m1,(LONG)r1*(LONG)2*c1*sizeof(double));
  462.         } 
  463.         else if (hd1->type==s_real || hd1->type==s_matrix)
  464.         {    make_complex(next_param(st));
  465.             mvconcat(st);
  466.             return;
  467.         }
  468.         else wrong_arg();
  469.     }    
  470.     else wrong_arg();
  471.     moveresult(st,result);
  472. }
  473.  
  474. void mhconcat (header *hd)
  475. {    header *st=hd,*hd1,*result;
  476.     double *m,*m1,*mr;
  477.     int r,c,r1,c1,i;
  478.     hd=getvalue(hd);
  479.     hd1=next_param(st);
  480.     if (!hd1) wrong_arg();
  481.     hd1=getvalue(hd1); if (error) return;
  482.     if (hd->type==s_string && hd1->type==s_string)
  483.     {    result=new_string(stringof(hd),
  484.             strlen(stringof(hd))+strlen(stringof(hd1))+1,"");
  485.         strcat(stringof(result),stringof(hd1));
  486.     }
  487.     else if (hd->type==s_real || hd->type==s_matrix)
  488.     {    if (hd1->type==s_real || hd1->type==s_matrix)
  489.         {    getmatrix(hd,&r,&c,&m);
  490.             getmatrix(hd1,&r1,&c1,&m1);
  491.             if (c!=0 && (r!=r1 || (LONG)c+c1>INT_MAX)) wrong_arg();
  492.             result=new_matrix(r1,c+c1,"");
  493.             mr=matrixof(result);
  494.             for (i=0; i<r1; i++)
  495.             {    if (c!=0) memmove((char *)mat(mr,c+c1,i,0),
  496.                     (char *)mat(m,c,i,0),c*sizeof(double));
  497.                 memmove((char *)mat(mr,c+c1,i,c),
  498.                     (char *)mat(m1,c1,i,0),c1*sizeof(double));
  499.             }
  500.         }
  501.         else if (hd1->type==s_complex || hd1->type==s_cmatrix)
  502.         {    make_complex(st);
  503.             mhconcat(st);
  504.             return;
  505.         }
  506.         else wrong_arg();
  507.     }
  508.     else if (hd->type==s_complex || hd->type==s_cmatrix)
  509.     {    if (hd1->type==s_complex || hd1->type==s_cmatrix)
  510.         {    getmatrix(hd,&r,&c,&m);
  511.             getmatrix(hd1,&r1,&c1,&m1);
  512.             if (c!=0 && (r!=r1 || (LONG)c+c1>INT_MAX)) wrong_arg();
  513.             result=new_cmatrix(r1,c+c1,"");
  514.             mr=matrixof(result);
  515.             for (i=0; i<r1; i++)
  516.             {    if (c!=0) memmove((char *)cmat(mr,c+c1,i,0),
  517.                     (char *)cmat(m,c,i,0),c*(LONG)2*sizeof(double));
  518.                 memmove((char *)cmat(mr,c+c1,i,c),
  519.                     (char *)cmat(m1,c1,i,0),c1*(LONG)2*sizeof(double));
  520.             }
  521.         }
  522.         else if (hd1->type==s_real || hd1->type==s_matrix)
  523.         {    make_complex(next_param(st)); if (error) return;
  524.             mhconcat(st);
  525.             return;
  526.         }
  527.         else wrong_arg();
  528.     }    
  529.     else wrong_arg();
  530.     moveresult(st,result);
  531. }
  532.  
  533. void cscalp (double *s, double *si, double *x, double *xi,
  534.     double *y, double *yi);
  535. void ccopy (double *y, double *x, double *xi);
  536.  
  537. void wmultiply (header *hd)
  538. /***** multiply
  539.     matrix multiplication for weakly nonzero matrices.
  540. *****/
  541. {    header *result,*st=hd,*hd1;
  542.     dims *d,*d1;
  543.     double *m,*m1,*m2,*mm1,*mm2,x,y,null=0.0;
  544.     int i,j,c,r,k;
  545.     hd=getvalue(hd); hd1=getvalue(nextof(st));
  546.     if (error) return;
  547.     if (hd->type==s_matrix && hd1->type==s_matrix)
  548.     {    d=dimsof(hd);
  549.         d1=dimsof(hd1);
  550.         if (d->c != d1->r)
  551.         {    error=8; output("Cannot multiply these!\n");
  552.             return;
  553.         }
  554.         r=d->r; c=d1->c;
  555.         result=new_matrix(r,c,"");
  556.         if (error) return;
  557.         m=matrixof(result);
  558.         m1=matrixof(hd);
  559.         m2=matrixof(hd1);
  560.         for (i=0; i<r; i++)
  561.             for (j=0; j<c; j++)
  562.             {    mm1=mat(m1,d->c,i,0); mm2=m2+j;
  563.                 x=0.0;
  564.                 for (k=0; k<d->c; k++)
  565.                 {    if ((*mm1!=0.0)&&(*mm2!=0.0)) x+=(*mm1)*(*mm2);
  566.                     mm1++; mm2+=d1->c;
  567.                 }
  568.                 *mat(m,c,i,j)=x;
  569.             }
  570.         moveresult(st,result);
  571.         return;
  572.     }
  573.     if (hd->type==s_matrix && hd1->type==s_cmatrix)
  574.     {    d=dimsof(hd);
  575.         d1=dimsof(hd1);
  576.         if (d->c != d1->r)
  577.         {    error=8; output("Cannot multiply these!\n");
  578.             return;
  579.         }
  580.         r=d->r; c=d1->c;
  581.         result=new_cmatrix(r,c,"");
  582.         if (error) return;
  583.         m=matrixof(result);
  584.         m1=matrixof(hd);
  585.         m2=matrixof(hd1);
  586.         for (i=0; i<r; i++)
  587.             for (j=0; j<c; j++)
  588.             {   mm1=mat(m1,d->c,i,0); mm2=m2+(LONG)2*j;
  589.                 x=0.0; y=0.0;
  590.                 for (k=0; k<d->c; k++)
  591.                 {    if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
  592.                             (*mm1!=0.0))
  593.                     cscalp(&x,&y,mm1,&null,mm2,mm2+1);
  594.                     mm1++; mm2+=2*d1->c;
  595.                 }
  596.                 ccopy(cmat(m,c,i,j),&x,&y);
  597.             }
  598.         moveresult(st,result);
  599.         return;
  600.     }
  601.     if (hd->type==s_cmatrix && hd1->type==s_matrix)
  602.     {    d=dimsof(hd);
  603.         d1=dimsof(hd1);
  604.         if (d->c != d1->r)
  605.         {    error=8; output("Cannot multiply these!\n");
  606.             return;
  607.         }
  608.         r=d->r; c=d1->c;
  609.         result=new_cmatrix(r,c,"");
  610.         if (error) return;
  611.         m=matrixof(result);
  612.         m1=matrixof(hd);
  613.         m2=matrixof(hd1);
  614.         for (i=0; i<r; i++)
  615.             for (j=0; j<c; j++)
  616.             {    mm1=cmat(m1,d->c,i,0); mm2=m2+j;
  617.                 x=0.0; y=0.0;
  618.                 for (k=0; k<d->c; k++)
  619.                 {    if ((*mm2!=0.0) &&
  620.                             (*mm1!=0.0 ||  *(mm1+1)!=0.0))
  621.                     cscalp(&x,&y,mm1,mm1+1,mm2,&null);
  622.                     mm1+=2; mm2+=d1->c;
  623.                 }
  624.                 ccopy(cmat(m,c,i,j),&x,&y);
  625.             }
  626.         moveresult(st,result);
  627.         return;
  628.     }
  629.     if (hd->type==s_cmatrix && hd1->type==s_cmatrix)
  630.     {    d=dimsof(hd);
  631.         d1=dimsof(hd1);
  632.         if (d->c != d1->r)
  633.         {    error=8; output("Cannot multiply these!\n");
  634.             return;
  635.         }
  636.         r=d->r; c=d1->c;
  637.         result=new_cmatrix(r,c,"");
  638.         if (error) return;
  639.         m=matrixof(result);
  640.         m1=matrixof(hd);
  641.         m2=matrixof(hd1);
  642.         for (i=0; i<r; i++)
  643.             for (j=0; j<c; j++)
  644.             {   mm1=cmat(m1,d->c,i,0); mm2=m2+(LONG)2*j;
  645.                 x=0.0; y=0.0;
  646.                 for (k=0; k<d->c; k++)
  647.                 {    if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
  648.                             (*mm1!=0.0 || *(mm1+1)!=0.0))
  649.                     cscalp(&x,&y,mm1,mm1+1,mm2,mm2+1);
  650.                     mm1+=2; mm2+=2*d1->c;
  651.                 }
  652.                 ccopy(cmat(m,c,i,j),&x,&y);
  653.             }
  654.         moveresult(st,result);
  655.         return;
  656.     }
  657.     else wrong_arg();
  658. }
  659.  
  660. void smultiply (header *hd)
  661. /***** multiply
  662.     matrix multiplication for weakly nonzero symmetric matrices.
  663. *****/
  664. {    header *result,*st=hd,*hd1;
  665.     dims *d,*d1;
  666.     double *m,*m1,*m2,*mm1,*mm2,x,y,null=0.0;
  667.     int i,j,c,r,k;
  668.     hd=getvalue(hd); hd1=getvalue(nextof(st));
  669.     if (error) return;
  670.     if (hd->type==s_matrix && hd1->type==s_matrix)
  671.     {    d=dimsof(hd);
  672.         d1=dimsof(hd1);
  673.         if (d->c != d1->r)
  674.         {    error=8; output("Cannot multiply these!\n");
  675.             return;
  676.         }
  677.         r=d->r; c=d1->c;
  678.         result=new_matrix(r,c,"");
  679.         if (error) return;
  680.         m=matrixof(result);
  681.         m1=matrixof(hd);
  682.         m2=matrixof(hd1);
  683.         for (i=0; i<r; i++)
  684.             for (j=i; j<c; j++)
  685.             {    mm1=mat(m1,d->c,i,0); mm2=m2+j;
  686.                 x=0.0;
  687.                 for (k=0; k<d->c; k++)
  688.                 {    if ((*mm1!=0.0)&&(*mm2!=0.0)) x+=(*mm1)*(*mm2);
  689.                     mm1++; mm2+=d1->c;
  690.                 }
  691.                 *mat(m,c,i,j)=x;
  692.                 *mat(m,c,j,i)=x;
  693.             }
  694.         moveresult(st,result);
  695.         return;
  696.     }
  697.     if (hd->type==s_matrix && hd1->type==s_cmatrix)
  698.     {    d=dimsof(hd);
  699.         d1=dimsof(hd1);
  700.         if (d->c != d1->r)
  701.         {    error=8; output("Cannot multiply these!\n");
  702.             return;
  703.         }
  704.         r=d->r; c=d1->c;
  705.         result=new_cmatrix(r,c,"");
  706.         if (error) return;
  707.         m=matrixof(result);
  708.         m1=matrixof(hd);
  709.         m2=matrixof(hd1);
  710.         for (i=0; i<r; i++)
  711.             for (j=i; j<c; j++)
  712.             {   mm1=mat(m1,d->c,i,0); mm2=m2+(LONG)2*j;
  713.                 x=0.0; y=0.0;
  714.                 for (k=0; k<d->c; k++)
  715.                 {    if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
  716.                             (*mm1!=0.0))
  717.                     cscalp(&x,&y,mm1,&null,mm2,mm2+1);
  718.                     mm1++; mm2+=2*d1->c;
  719.                 }
  720.                 ccopy(cmat(m,c,i,j),&x,&y); y=-y;
  721.                 ccopy(cmat(m,c,j,i),&x,&y);
  722.             }
  723.         moveresult(st,result);
  724.         return;
  725.     }
  726.     if (hd->type==s_cmatrix && hd1->type==s_matrix)
  727.     {    d=dimsof(hd);
  728.         d1=dimsof(hd1);
  729.         if (d->c != d1->r)
  730.         {    error=8; output("Cannot multiply these!\n");
  731.             return;
  732.         }
  733.         r=d->r; c=d1->c;
  734.         result=new_cmatrix(r,c,"");
  735.         if (error) return;
  736.         m=matrixof(result);
  737.         m1=matrixof(hd);
  738.         m2=matrixof(hd1);
  739.         for (i=0; i<r; i++)
  740.             for (j=i; j<c; j++)
  741.             {    mm1=cmat(m1,d->c,i,0); mm2=m2+j;
  742.                 x=0.0; y=0.0;
  743.                 for (k=0; k<d->c; k++)
  744.                 {    if ((*mm2!=0.0) &&
  745.                             (*mm1!=0.0 ||  *(mm1+1)!=0.0))
  746.                     cscalp(&x,&y,mm1,mm1+1,mm2,&null);
  747.                     mm1+=2; mm2+=d1->c;
  748.                 }
  749.                 ccopy(cmat(m,c,i,j),&x,&y); y=-y;
  750.                 ccopy(cmat(m,c,j,i),&x,&y);
  751.             }
  752.         moveresult(st,result);
  753.         return;
  754.     }
  755.     if (hd->type==s_cmatrix && hd1->type==s_cmatrix)
  756.     {    d=dimsof(hd);
  757.         d1=dimsof(hd1);
  758.         if (d->c != d1->r)
  759.         {    error=8; output("Cannot multiply these!\n");
  760.             return;
  761.         }
  762.         r=d->r; c=d1->c;
  763.         result=new_cmatrix(r,c,"");
  764.         if (error) return;
  765.         m=matrixof(result);
  766.         m1=matrixof(hd);
  767.         m2=matrixof(hd1);
  768.         for (i=0; i<r; i++)
  769.             for (j=i; j<c; j++)
  770.             {   mm1=cmat(m1,d->c,i,0); mm2=m2+(LONG)2*j;
  771.                 x=0.0; y=0.0;
  772.                 for (k=0; k<d->c; k++)
  773.                 {    if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
  774.                             (*mm1!=0.0 || *(mm1+1)!=0.0))
  775.                     cscalp(&x,&y,mm1,mm1+1,mm2,mm2+1);
  776.                     mm1+=2; mm2+=2*d1->c;
  777.                 }
  778.                 ccopy(cmat(m,c,i,j),&x,&y); 
  779.                 y=-y;
  780.                 ccopy(cmat(m,c,j,i),&x,&y);
  781.             }
  782.         moveresult(st,result);
  783.         return;
  784.     }
  785.     else wrong_arg();
  786. }
  787.  
  788.