1 /*****************************************************************************
2  *
3  *  Elmer, A Finite Element Software for Multiphysical Problems
4  *
5  *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6  *
7  *  This program is free software; you can redistribute it and/or
8  *  modify it under the terms of the GNU General Public License
9  *  as published by the Free Software Foundation; either version 2
10  *  of the License, or (at your option) any later version.
11  *
12  *  This program is distributed in the hope that it will be useful,
13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *  GNU General Public License for more details.
16  *
17  *  You should have received a copy of the GNU General Public License
18  *  along with this program (in file fem/GPL-2); if not, write to the
19  *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
20  *  Boston, MA 02110-1301, USA.
21  *
22  *****************************************************************************/
23 
24 /*******************************************************************************
25  *
26  *     Evaluate expression threes in format parsed by parser.c.
27  *
28  *******************************************************************************
29  *
30  *                     Author:       Juha Ruokolainen
31  *
32  *                    Address: CSC - IT Center for Science Ltd.
33  *                                Keilaranta 14, P.O. BOX 405
34  *                                  02101 Espoo, Finland
35  *                                  Tel. +358 0 457 2723
36  *                                Telefax: +358 0 457 2302
37  *                              EMail: Juha.Ruokolainen@csc.fi
38  *
39  *                       Date: 30 May 1996
40  *
41  *                Modified by:
42  *
43  *       Date of modification:
44  *
45  ******************************************************************************/
46 /***********************************************************************
47 |
48 |  EVAL.C - Last Edited 6. 8. 1988
49 |
50 ***********************************************************************/
51 #include "../config.h"
52 /*======================================================================
53 |Syntax of the manual pages:
54 |
55 |FUNCTION NAME(...) params ...
56 |
57 $  usage of the function and type of the parameters
58 ?  explane the effects of the function
59 =  return value and the type of value if not of type int
60 @  globals effected directly by this routine
61 !  current known bugs or limitations
62 &  functions called by this function
63 ~  these functions may interest you as an alternative function or
64 |  because they control this function somehow
65 ^=====================================================================*/
66 
67 
68 /*
69  * $Id: eval.c,v 1.5 2006/02/02 06:51:16 jpr Exp $
70  *
71  * $Log: eval.c,v $
72  * Revision 1.4  2005/08/25 13:44:22  vierinen
73  * windoze stuff
74  *
75  * Revision 1.3  2005/05/27 12:26:20  vierinen
76  * changed header install location
77  *
78  * Revision 1.2  2005/05/26 12:34:53  vierinen
79  * windows stuff
80  *
81  * Revision 1.1.1.1  2005/04/14 13:29:14  vierinen
82  * initial matc automake package
83  *
84  * Revision 1.2  1998/08/01 12:34:35  jpr
85  *
86  * Added Id, started Log.
87  *
88  *
89  */
90 
91 #include "elmer/matc.h"
92 
evaltree(root)93 VARIABLE *evaltree(root) TREE *root;
94 /*======================================================================
95 ?  Evaluate the equation tree (real hard). The tree might be a list of
96 |  trees. If there is more than one tree in the list, the return value
97 |  will be a single column vector, the result from each tree is
98 |  concatenated to this vector. If there's only one tree, return
99 |  value is tree's evaluated result.
100 |
101 &  strlen(), fopen(), fclose(), var_temp_new(), var_temp_copy(),
102 |  var_delete_temp(), com_check(), var_check(), fnc_check(),
103 |  fnc_exec(), com_pointw(), com_source()
104 |
105 ^=====================================================================*/
106 {
107 
108   VARIABLE *first,     /* used to hold start address for    */
109                        /* temporary list of VARIABLES, the  */
110                        /* evaluated trees results           */
111             *subs,     /* indexes for tree elements         */
112             *par,      /* parameters for tree elements      */
113             *res,      /* result is returned in this var.   */
114                        /* probably used as a temporary too  */
115             *leftptr,  /* result from left leaf             */
116             *rightptr, /* result from right leaf            */
117             *tmp,      /* can't manage without these        */
118             *tmp1;
119 
120   FUNCTION *fnc;       /* needed, if someone's using defined functions */
121 
122   COMMAND *com;        /* needed, maybe */
123 
124   MATRIX *opres;       /* don't know really */
125 
126   TREE *argptr;        /* temporary */
127 
128   int i,               /* there's always reason for these         */
129       dim,             /* the final dimension of the return value */
130       argcount;        /* parameter count                         */
131 
132   char *resbeg;        /* for memcpying */
133 
134   FILE *fp;            /* used to check if a file exists */
135 
136 
137   if (root == NULL) return NULL;
138 
139   dim = 0;  first = NULL;
140 
141 /*-------------------------------------------------------------------*/
142 
143   /*
144       while there's trees in the list.
145   */
146   while(root)
147   {
148 
149     subs = par = tmp = NULL;
150 
151     /*
152         check if the result will be indexed
153     */
154     argptr = SUBS(root);
155     if (argptr)
156     {
157       subs = tmp = evaltree(argptr);
158       argptr = NEXT(argptr);
159       while(argptr)
160       {
161         NEXT(tmp) = evaltree(argptr);
162         argptr = NEXT(argptr); tmp = NEXT(tmp);
163       }
164     }
165 
166     switch(ETYPE(root))
167     {
168 
169     /******************************************************
170               some kind of existing identifier.
171     *******************************************************/
172     case ETYPE_NAME:
173       /*
174           is there parameters for this one ?
175           and how many ?
176       */
177       argptr = ARGS(root);
178       argcount = 0;
179       if (argptr)
180       {
181         par = tmp = evaltree(argptr);
182         argptr = NEXT(argptr);
183         argcount++;
184         while(argptr)
185         {
186           argcount++;
187           NEXT(tmp) = evaltree(argptr);
188           argptr = NEXT(argptr); tmp = NEXT(tmp);
189         }
190       }
191 
192 
193       /*
194           one of the builtin commands (sin, ones, help, ...)
195       */
196       if ((com = com_check(SDATA(root))) != (COMMAND *)NULL)
197       {
198 
199         if (argcount < com->minp || argcount > com->maxp)
200         {
201           if (com->minp == com->maxp)
202           {
203             error( "Builtin function [%s] requires %d argument(s).\n", SDATA(root), com->minp);
204           }
205           else
206           {
207              error("Builtin function [%s] takes from %d to %d argument(s).\n",
208                        SDATA(root), com->minp, com->maxp);
209           }
210         }
211 
212 
213         if (com->flags & CMDFLAG_PW)
214         {
215           tmp = com_pointw((double (*)())com->sub, par);
216         }
217         else
218         {
219           tmp = (*com->sub)(par);
220         }
221       }
222 
223 
224       /*
225           a variable name
226       */
227       else if ((tmp1 = var_check(SDATA(root))) != (VARIABLE *)NULL)
228       {
229         tmp = (VARIABLE *)ALLOCMEM(VARIABLESIZE);
230         tmp ->this = tmp1->this;
231         REFCNT(tmp)++;
232         if (par != NULL)
233         {
234           subs = par; par = NULL;
235         }
236       }
237 
238       /*
239           user defined function
240       */
241       else if ((fnc = fnc_check(SDATA(root))) != (FUNCTION *)NULL)
242       {
243         tmp = fnc_exec(fnc, par); par = NULL;
244       }
245 
246 
247       /*
248            maybe a file name ?
249       */
250       else if ((fp = fopen(SDATA(root),"r")) != (FILE *)NULL)
251       {
252         fclose(fp);
253         tmp = var_temp_new(TYPE_STRING, 1, strlen(SDATA(root)));
254         for(i = 0; i < strlen(SDATA(root)); i++)
255         {
256           M(tmp, 0, i) = SDATA(root)[i];
257         }
258         (*com_source)(tmp);
259         var_delete_temp(tmp);
260         tmp = NULL;
261       }
262 
263       /*
264           troubles!
265       */
266       else
267       {
268         error("Undeclared identifier: [%s].\n", SDATA(root));
269       }
270 
271       break;
272 
273 
274     /******************************************************
275                   single string constant
276     *******************************************************/
277     case ETYPE_STRING:
278       tmp = var_temp_new(TYPE_STRING, 1, strlen(SDATA(root)));
279       for(i = 0; i < strlen(SDATA(root)); i++)
280       {
281         M(tmp,0,i) = SDATA(root)[i];
282       }
283       break;
284 
285     /******************************************************
286                   single numeric constant
287     *******************************************************/
288     case ETYPE_NUMBER:
289       tmp = var_temp_new(TYPE_DOUBLE, 1, 1);
290       M(tmp,0,0) = DDATA(root);
291       break;
292 
293     /******************************************************
294                     constant matrix
295     *******************************************************/
296     case ETYPE_CONST:
297       tmp = (VARIABLE *)ALLOCMEM(VARIABLESIZE);
298       tmp->this = CDATA(root)->this;
299       REFCNT(tmp)++;
300       break;
301 
302     /******************************************************
303                            huh ?
304     *******************************************************/
305     case ETYPE_EQUAT:
306       tmp = evaltree(LEFT(root));
307       break;
308 
309     /******************************************************
310          left oper [right]
311          oper = divide, multiply, transpose, power,...
312     *******************************************************/
313     case ETYPE_OPER:
314       /*
315        *   evaluate both leafs.
316        */
317       leftptr = rightptr = NULL; opres = NULL;
318 
319       leftptr  = evaltree(LEFT(root));
320       rightptr = evaltree(RIGHT(root));
321 
322       if (leftptr && rightptr)
323         opres = (*VDATA(root))(leftptr->this, rightptr->this);
324       else if (leftptr)
325         opres = (*VDATA(root))(leftptr->this, NULL);
326       else if (rightptr)
327         opres = (*VDATA(root))(rightptr->this, NULL);
328 
329       var_delete_temp(leftptr);
330       var_delete_temp(rightptr);
331 
332       if (opres != NULL)
333       {
334         tmp = (VARIABLE *)ALLOCMEM(VARIABLESIZE);
335         tmp->this = opres;
336         REFCNT(tmp) = 1;
337       }
338 
339       break;
340     }
341 
342     /*
343        if NULL result, don't really know what to do, so we try quitting.
344     */
345     /*
346     if (tmp == (VARIABLE *)NULL)
347     {
348       if (subs) var_delete_temp(subs);
349       if (par) var_delete_temp(par);
350       if (first) var_delete_temp(first);
351       return (VARIABLE *)NULL;
352     }
353     */
354 
355     /******************************************************
356        deleteting temporaries, preparing for next loop
357     *******************************************************/
358     if (subs != NULL)
359     {
360       if (tmp != NULL)
361       {
362         tmp1 = tmp;
363         NEXT(tmp1) = subs;
364         tmp = com_el(tmp1);
365         var_delete_temp(tmp1);
366       }
367       else
368       {
369         var_delete_temp(subs);
370       }
371       tmp1 = subs = NULL;
372     }
373 
374     if (first == NULL)
375     {
376       first = res = tmp;
377     }
378     else if (tmp != NULL)
379     {
380       NEXT(res) = tmp; res = NEXT(res);
381     }
382 
383     if (subs) var_delete_temp(subs);
384     if (par) var_delete_temp(par);
385 
386     if (tmp != NULL) dim += NROW(tmp) * NCOL(tmp);/* count size of the result */
387 
388     root = LINK(root);
389   }
390  /*-------------------------------------------------------------------*/
391 
392   /*
393       there was only one tree, so return  it's result.
394   */
395   if (tmp == first) return first;
396 
397   /*
398       it was a list of trees, concatenate the results
399   */
400   res = var_temp_new(TYPE(first), 1, dim);
401 
402   resbeg = (char *)MATR(res);
403   for(tmp = first; tmp; tmp = NEXT(tmp))
404   {
405     memcpy(resbeg, (char *)MATR(tmp), MATSIZE(tmp));
406     resbeg += MATSIZE(tmp);
407   }
408   /*
409       and delete rest of the temporaries.
410   */
411   var_delete_temp(first);
412 
413   return res;
414 }
415 
evalclause(root)416 VARIABLE *evalclause(root) CLAUSE *root;
417 /*======================================================================
418 ?  Evaluate the operations list. The list contains equations trees
419 |  (actually lists of trees :-).
420 |
421 &  ALLOCMEM, FREEMEM, STRCOPY, var_temp_new(), var_delete_temp(),
422 &  var_check(), com_check(), fnc_check(), evaltree(),
423 ^=====================================================================*/
424 {
425   VARIABLE *ptr = NULL,
426            *res, *par, *tmp;      /* used for something magic */
427 
428   TREE *argptr;                   /* pointer to parameters    */
429 
430   double *d;
431 
432   int i;
433 
434 /*-------------------------------------------------------------------*/
435 
436   while(root)
437   {
438     if (root->data == endsym) return ptr;
439 
440     switch(root->data)
441     {
442       /************************************************************
443                        System call
444       ************************************************************/
445       case systemcall:
446       {
447 #if defined(WIN32) || defined(MINGW32)
448            FILE *fp = _popen( SDATA(root->this), "r" );
449 #else
450            FILE *fp = popen( SDATA(root->this), "r" );
451 #endif
452            static char s[101];
453 
454            if ( !fp ) error( "systemcall: open failure: [%s].\n", SDATA(root->this) );
455 
456            while( fgets( s, 120, fp ) ) PrintOut( s );
457 
458 #if defined(WIN32) || defined(MINGW32)
459            _pclose( fp );
460 #else
461            pclose( fp );
462 #endif
463       }
464       break;
465       /************************************************************
466                        Function definition
467       ************************************************************/
468       case funcsym:
469       {
470 
471         FUNCTION *fnc;                    /* pointer to created function */
472         TREE *tptr;                       /*  function parameters        */
473         char *name = SDATA(root->this);   /* function name */
474         int i, n, argcount;               /* well... */
475 
476         /*
477             check for name conflicts
478         */
479         if (var_check(name) || com_check(name))
480         {
481           error( "Function not created [%s], identifier in use.\n",name);
482         }
483 
484         /*
485          *   allocate mem for FUNCTION structure and add it
486          *   to the the FUNCTIONS list
487          */
488         if (fnc = fnc_check(name)) fnc_free_entry(fnc);
489         fnc = (FUNCTION *)ALLOCMEM(sizeof(FUNCTION));
490         NAME(fnc) = STRCOPY(name);
491         lst_add(FUNCTIONS, (LIST *)fnc);
492 
493         /*
494          *   copy parameter names to the structure.
495          */
496         argcount = 0;
497         for(tptr = ARGS(root->this); tptr; tptr = NEXT(tptr)) argcount++;
498         if (argcount > 0)
499         {
500           fnc -> parnames = (char **)ALLOCMEM(argcount * sizeof(char *));
501           for(i = 0, tptr = ARGS(root->this); tptr; tptr = NEXT(tptr))
502             fnc -> parnames[i++] = STRCOPY(SDATA(tptr));
503         }
504         fnc -> parcount = argcount;
505 
506         /*
507          *   copy help text if any to the structure.
508          */
509         argcount = n = 0;
510         for( tptr = SUBS(root->this); tptr; tptr = NEXT(tptr))
511         {
512             if ( SDATA(tptr) )
513             {
514                 argcount++;
515                 n += strlen( SDATA(tptr) );
516             }
517         }
518         if ( argcount > 0 && n > 0)
519         {
520             fnc->help = (char *)ALLOCMEM( (n+argcount+1)*sizeof(char) );
521             for( tptr = SUBS(root->this); tptr; tptr = NEXT(tptr) )
522             {
523                 if ( SDATA(tptr) )
524                 {
525                     strcat( fnc->help, SDATA(tptr) );
526                     strcat( fnc->help, "\n" );
527                 }
528             }
529         }
530 
531         /*
532          *   copy imported variable names to the structure.
533          */
534         argcount = 0;
535         for(tptr = LEFT(root->this); tptr; tptr = NEXT(tptr)) argcount++;
536         if (argcount > 0)
537         {
538           fnc -> imports = (char **)ALLOCMEM((argcount+1) * sizeof(char *));
539           for(i = 0, tptr = LEFT(root->this); tptr; tptr = NEXT(tptr))
540             fnc -> imports[i++] = STRCOPY(SDATA(tptr));
541           fnc -> imports[i] = NULL;
542         }
543         else
544           fnc -> imports = NULL;
545 
546         /*
547          *   copy exported variable names to the structure.
548          */
549         argcount = 0;
550         for(tptr = RIGHT(root->this); tptr; tptr = NEXT(tptr)) argcount++;
551         if (argcount > 0)
552         {
553           fnc -> exports = (char **)ALLOCMEM((argcount+1) * sizeof(char *));
554           for(i = 0, tptr = RIGHT(root->this); tptr; tptr = NEXT(tptr))
555             fnc -> exports[i++] = STRCOPY(SDATA(tptr));
556           fnc -> exports[i] = NULL;
557         }
558         else
559           fnc -> exports = NULL;
560 
561         /*
562             and finally the function body
563         */
564         fnc -> body = LINK(root); LINK(root) = NULL;
565 
566 /*        PrintOut( "Function defined: [%s].\n", name); */
567 
568         return NULL;
569       }
570 
571 
572       /***************************************************************
573                              statement
574       ****************************************************************/
575       case assignsym:
576       {
577 
578         int  iflg = FALSE,          /* is the result indexed   */
579              pflg = TRUE;           /* should it be printed to */
580                                     /* output stream           */
581 
582         char *r = "ans";            /* resulted VARIABLE name  */
583 
584         par = NULL;
585 
586         /*
587          *   there is an explicit assigment  (ie. x = x + 1, not x + 1)
588          */
589         if (root->this)
590         {
591           /*
592            *  VARIABLE name
593            */
594           r = SDATA( root->this );
595 
596           /*
597            *  check for name conflicts
598            */
599           if ( fnc_check(r) || com_check(r) || lst_find(CONSTANTS, r) )
600           {
601               error( "VARIABLE not created [%s], identifier in use.\n", r );
602           }
603 
604           /*
605            *   is it indexed ?
606            */
607           pflg = FALSE;
608           argptr = ARGS(root->this);
609           if (argptr)
610           {
611             iflg = TRUE;
612             if ((par = tmp = evaltree(argptr)) != NULL)
613             {
614               argptr = NEXT(argptr);
615               while(argptr)
616               {
617                 if ((NEXT(tmp) = evaltree(argptr)) == NULL) break;
618                 argptr = NEXT(argptr);
619                 tmp = NEXT(tmp);
620               }
621             }
622           }
623         }
624 
625         /*
626          *  evaluate the right side of the statement
627          *  and put the result where it belongs
628          */
629         ptr = evaltree( LINK(root)->this );
630         ptr = put_result( ptr, r, par, iflg, pflg );
631 
632         if ( par ) var_delete_temp( par );
633         root = LINK(root);
634 
635       break;
636       }
637 
638       /***************************************************************
639                              if statement
640       ****************************************************************/
641       case ifsym:
642 
643         if ((res = evaltree(root->this)) != NULL)
644         {
645           for(i = 0, d=MATR(res); i < NROW(res)*NCOL(res); i++)
646             if (*d++ == 0) break;
647 
648           /*
649            *   condition false
650            */
651           if (*--d == 0)
652           {
653             root = root->jmp;
654             if (root->data == elsesym)
655             {
656               ptr = evalclause(LINK(root));
657               root = root -> jmp;
658             }
659           }
660 
661           /*
662            *   condition true
663            */
664           else
665           {
666             ptr = evalclause(LINK(root));
667             root = root->jmp;
668             if (root->data == elsesym)
669             {
670               root = root->jmp;
671             }
672           }
673           var_delete_temp(res);
674         }
675         else
676         {
677           root = root -> jmp;
678           if (root->data == elsesym)
679           {
680             root = root->jmp;
681           }
682         }
683       break;
684 
685       /***************************************************************
686                              while statement
687       ****************************************************************/
688       case whilesym:
689 
690         while(TRUE)
691         {
692           if ((res = evaltree(root->this)) == NULL) break;
693 
694           for(i=0, d=MATR(res); i < NROW(res)*NCOL(res); i++)
695             if (*d++ == 0) break;
696 
697           /*
698               condition true, go for another loop
699           */
700           if (*--d != 0)
701           {
702             ptr = evalclause(LINK(root));
703             var_delete_temp(res);
704           }
705 
706           /*
707               condition false, done with while-loop
708           */
709           else
710           {
711             var_delete_temp(res);
712             break;
713           }
714         }
715         root = root->jmp;
716       break;
717 
718       /***************************************************************
719                              for statement
720       ****************************************************************/
721       case forsym:
722       {
723         VARIABLE *var;
724         char *r;
725 
726         /*
727          *  VARIABLE name
728          */
729         r = SDATA(root->this);
730 
731         /*
732          *  check for name conflicts
733          */
734         if (fnc_check(r) || com_check(r) || lst_find(CONSTANTS, r))
735         {
736           error( "VARIABLE not created [%s], identifier in use.\n ", r);
737         }
738 
739         if ((res = evaltree(LINK(root->this))) != NULL)
740         {
741           var_delete(r);
742           var = var_new(r,TYPE(res),1,1);
743 
744           d = MATR(res);
745           for(i = 0; i < NCOL(res)*NROW(res); i++)
746           {
747             *MATR(var) = *d++;
748             ptr = evalclause(LINK(root));
749           }
750 
751           var_delete_temp(res);
752         }
753         root = root->jmp;
754       break;
755       }
756     }
757     root = LINK(root);
758   }
759   return ptr;
760 }
761 
put_values(ptr,resname,par)762 VARIABLE *put_values(ptr, resname, par) VARIABLE *ptr, *par; char *resname;
763 /*======================================================================
764 ?  extract values from ptr, indexes in par, and put them to res
765 ^=====================================================================*/
766 {
767   static double defind = 0.0;
768 
769   double *ind1,
770          *ind2,
771          *dtmp;
772 
773   int i, j, k,
774       rows, cols,
775       size1, size2,
776       imax1, imax2,
777       ind;
778 
779   VARIABLE *res;
780 
781   res = (VARIABLE *)lst_find(VARIABLES, resname);
782 
783   if (NEXT(par) == NULL)
784   {
785     if (
786       res != NULL && NROW(par) == NROW(res) && NCOL(par) == NCOL(res)
787                   && !(NROW(res) == 1 && NCOL(res) == 1)
788       )
789     {
790       int logical = TRUE,
791           csize   = 0;
792 
793       dtmp = MATR(par);
794       for(i = 0; i < NROW(par)*NCOL(par); i++)
795         if (dtmp[i] != 0 && dtmp[i] != 1)
796         {
797           logical = FALSE;
798           break;
799         }
800 
801       if (logical)
802       {
803         imax1 = NROW(ptr) * NCOL(ptr);
804         dtmp  = MATR(ptr);
805         for(i = 0, k = 0; i < NROW(res); i++)
806           for(j = 0,csize=0; j < NCOL(res); j++)
807           {
808             while(M(par,i,j)==1 && j+csize<NCOL(res) && k+csize<imax1) csize++;
809             if (csize > 0)
810             {
811               memcpy(&M(res,i,j),&dtmp[k],csize*sizeof(double));
812               j += csize-1;
813               k += csize;
814               csize = 0;
815               if (k >= imax1) k = 0;
816             }
817           }
818          var_delete_temp(ptr);
819          return res;
820       }
821       else
822       {
823         ind1  = &defind;
824         size1 = 1;
825         ind2  = MATR(par);
826         size2 = NCOL(par);
827       }
828     }
829     else
830     {
831       ind1  = &defind;
832       size1 = 1;
833       ind2  = MATR(par);
834       size2 = NCOL(par);
835     }
836   }
837   else
838   {
839     ind1  = MATR(par);
840     size1 = NCOL(par);
841     ind2  = MATR(NEXT(par));
842     size2 = NCOL(NEXT(par));
843   }
844 
845   imax1 = (int)ind1[0];
846   for(i = 1; i < size1; i++)
847   {
848     imax1 = max(imax1, (int)ind1[i]);
849   }
850 
851   imax2 = (int)ind2[0];
852   for(i = 1; i < size2; i++)
853   {
854     imax2 = max(imax2, (int)ind2[i]);
855   }
856 
857   if (res == NULL)
858   {
859     res = var_new(resname, TYPE(ptr), imax1+1, imax2+1);
860   }
861   else if (NROW(res) <= imax1 || NCOL(res) <= imax2)
862   {
863     int ir = NROW(res), jc = NCOL(res);
864     MATRIX *t;
865 
866     imax1 = max(ir, imax1 + 1);
867     imax2 = max(jc, imax2 + 1);
868 
869     t = mat_new(TYPE(res), imax1, imax2);
870     dtmp = t->data;
871 
872     for(i = 0; i < ir; i++)
873     {
874       memcpy(&dtmp[i*imax2],&M(res,i,0),jc*sizeof(double));
875     }
876 
877     if (--REFCNT(res) == 0)
878     {
879       mat_free(res->this);
880     }
881     res->this = t;
882     REFCNT(res) = 1;
883   }
884   else if (REFCNT(res) > 1)
885   {
886     --REFCNT(res);
887     res->this = mat_copy(res->this);
888   }
889 
890   imax1 = NROW(ptr) * NCOL(ptr);
891   dtmp = MATR(ptr);
892   for(i = 0, k = 0; i < size1; i++)
893   {
894     ind = (int)ind1[i];
895     for(j = 0; j < size2; j++)
896     {
897       memcpy(&M(res,ind,(int)ind2[j]),&dtmp[k++],sizeof(double));
898       if (k >= imax1) k = 0;
899     }
900   }
901 
902   var_delete_temp(ptr);
903 
904   return res;
905 }
906 
put_result(ptr,resname,par,indexflag,printflag)907 VARIABLE *put_result(ptr, resname, par, indexflag, printflag)
908 /*======================================================================
909 ?  copy VARIABLE from one place to another
910 |  and conditionally print the result
911 ^=====================================================================*/
912    int indexflag, printflag;
913    VARIABLE *ptr, *par;
914    char *resname;
915 {
916    VARIABLE *res = NULL;
917 
918   var_delete( "ans" );
919   if (indexflag && par)
920   {
921     res = put_values(ptr, resname, par);
922   }
923   else
924   {
925     res = var_rename( ptr, resname );
926   }
927 
928   if ( res ) res->changed = 1;
929   if (printflag) var_print(res);
930 
931   return res;
932 }
933