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 library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * This library 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 GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with this library (in file ../LGPL-2.1); if not, write
19  * to the Free Software Foundation, Inc., 51 Franklin Street,
20  * Fifth Floor, 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 #pragma omp threadprivate (s)
454 
455            if ( !fp ) error( "systemcall: open failure: [%s].\n", SDATA(root->this) );
456 
457            while( fgets( s, 120, fp ) ) PrintOut( s );
458 
459 #if defined(WIN32) || defined(MINGW32)
460            _pclose( fp );
461 #else
462            pclose( fp );
463 #endif
464       }
465       break;
466       /************************************************************
467                        Function definition
468       ************************************************************/
469       case funcsym:
470       {
471 
472         FUNCTION *fnc;                    /* pointer to created function */
473         TREE *tptr;                       /*  function parameters        */
474         char *name = SDATA(root->this);   /* function name */
475         int i, n, argcount;               /* well... */
476 
477         /*
478             check for name conflicts
479         */
480         if (var_check(name) || com_check(name))
481         {
482           error( "Function not created [%s], identifier in use.\n",name);
483         }
484 
485         /*
486          *   allocate mem for FUNCTION structure and add it
487          *   to the the FUNCTIONS list
488          */
489         if (fnc = fnc_check(name)) fnc_free_entry(fnc);
490         fnc = (FUNCTION *)ALLOCMEM(sizeof(FUNCTION));
491         NAME(fnc) = STRCOPY(name);
492         lst_add(FUNCTIONS, (LIST *)fnc);
493 
494         /*
495          *   copy parameter names to the structure.
496          */
497         argcount = 0;
498         for(tptr = ARGS(root->this); tptr; tptr = NEXT(tptr)) argcount++;
499         if (argcount > 0)
500         {
501           fnc -> parnames = (char **)ALLOCMEM(argcount * sizeof(char *));
502           for(i = 0, tptr = ARGS(root->this); tptr; tptr = NEXT(tptr))
503             fnc -> parnames[i++] = STRCOPY(SDATA(tptr));
504         }
505         fnc -> parcount = argcount;
506 
507         /*
508          *   copy help text if any to the structure.
509          */
510         argcount = n = 0;
511         for( tptr = SUBS(root->this); tptr; tptr = NEXT(tptr))
512         {
513             if ( SDATA(tptr) )
514             {
515                 argcount++;
516                 n += strlen( SDATA(tptr) );
517             }
518         }
519         if ( argcount > 0 && n > 0)
520         {
521             fnc->help = (char *)ALLOCMEM( (n+argcount+1)*sizeof(char) );
522             for( tptr = SUBS(root->this); tptr; tptr = NEXT(tptr) )
523             {
524                 if ( SDATA(tptr) )
525                 {
526                     strcat( fnc->help, SDATA(tptr) );
527                     strcat( fnc->help, "\n" );
528                 }
529             }
530         }
531 
532         /*
533          *   copy imported variable names to the structure.
534          */
535         argcount = 0;
536         for(tptr = LEFT(root->this); tptr; tptr = NEXT(tptr)) argcount++;
537         if (argcount > 0)
538         {
539           fnc -> imports = (char **)ALLOCMEM((argcount+1) * sizeof(char *));
540           for(i = 0, tptr = LEFT(root->this); tptr; tptr = NEXT(tptr))
541             fnc -> imports[i++] = STRCOPY(SDATA(tptr));
542           fnc -> imports[i] = NULL;
543         }
544         else
545           fnc -> imports = NULL;
546 
547         /*
548          *   copy exported variable names to the structure.
549          */
550         argcount = 0;
551         for(tptr = RIGHT(root->this); tptr; tptr = NEXT(tptr)) argcount++;
552         if (argcount > 0)
553         {
554           fnc -> exports = (char **)ALLOCMEM((argcount+1) * sizeof(char *));
555           for(i = 0, tptr = RIGHT(root->this); tptr; tptr = NEXT(tptr))
556             fnc -> exports[i++] = STRCOPY(SDATA(tptr));
557           fnc -> exports[i] = NULL;
558         }
559         else
560           fnc -> exports = NULL;
561 
562         /*
563             and finally the function body
564         */
565         fnc -> body = LINK(root); LINK(root) = NULL;
566 
567 /*        PrintOut( "Function defined: [%s].\n", name); */
568 
569         return NULL;
570       }
571 
572 
573       /***************************************************************
574                              statement
575       ****************************************************************/
576       case assignsym:
577       {
578 
579         int  iflg = FALSE,          /* is the result indexed   */
580              pflg = TRUE;           /* should it be printed to */
581                                     /* output stream           */
582 
583         char *r = "ans";            /* resulted VARIABLE name  */
584 
585         par = NULL;
586 
587         /*
588          *   there is an explicit assigment  (ie. x = x + 1, not x + 1)
589          */
590         if (root->this)
591         {
592           /*
593            *  VARIABLE name
594            */
595           r = SDATA( root->this );
596 
597           /*
598            *  check for name conflicts
599            */
600           if ( fnc_check(r) || com_check(r) || lst_find(CONSTANTS, r) )
601           {
602               error( "VARIABLE not created [%s], identifier in use.\n", r );
603           }
604 
605           /*
606            *   is it indexed ?
607            */
608           pflg = FALSE;
609           argptr = ARGS(root->this);
610           if (argptr)
611           {
612             iflg = TRUE;
613             if ((par = tmp = evaltree(argptr)) != NULL)
614             {
615               argptr = NEXT(argptr);
616               while(argptr)
617               {
618                 if ((NEXT(tmp) = evaltree(argptr)) == NULL) break;
619                 argptr = NEXT(argptr);
620                 tmp = NEXT(tmp);
621               }
622             }
623           }
624         }
625 
626         /*
627          *  evaluate the right side of the statement
628          *  and put the result where it belongs
629          */
630         ptr = evaltree( LINK(root)->this );
631         ptr = put_result( ptr, r, par, iflg, pflg );
632 
633         if ( par ) var_delete_temp( par );
634         root = LINK(root);
635 
636       break;
637       }
638 
639       /***************************************************************
640                              if statement
641       ****************************************************************/
642       case ifsym:
643 
644         if ((res = evaltree(root->this)) != NULL)
645         {
646           for(i = 0, d=MATR(res); i < NROW(res)*NCOL(res); i++)
647             if (*d++ == 0) break;
648 
649           /*
650            *   condition false
651            */
652           if (*--d == 0)
653           {
654             root = root->jmp;
655             if (root->data == elsesym)
656             {
657               ptr = evalclause(LINK(root));
658               root = root -> jmp;
659             }
660           }
661 
662           /*
663            *   condition true
664            */
665           else
666           {
667             ptr = evalclause(LINK(root));
668             root = root->jmp;
669             if (root->data == elsesym)
670             {
671               root = root->jmp;
672             }
673           }
674           var_delete_temp(res);
675         }
676         else
677         {
678           root = root -> jmp;
679           if (root->data == elsesym)
680           {
681             root = root->jmp;
682           }
683         }
684       break;
685 
686       /***************************************************************
687                              while statement
688       ****************************************************************/
689       case whilesym:
690 
691         while(TRUE)
692         {
693           if ((res = evaltree(root->this)) == NULL) break;
694 
695           for(i=0, d=MATR(res); i < NROW(res)*NCOL(res); i++)
696             if (*d++ == 0) break;
697 
698           /*
699               condition true, go for another loop
700           */
701           if (*--d != 0)
702           {
703             ptr = evalclause(LINK(root));
704             var_delete_temp(res);
705           }
706 
707           /*
708               condition false, done with while-loop
709           */
710           else
711           {
712             var_delete_temp(res);
713             break;
714           }
715         }
716         root = root->jmp;
717       break;
718 
719       /***************************************************************
720                              for statement
721       ****************************************************************/
722       case forsym:
723       {
724         VARIABLE *var;
725         char *r;
726 
727         /*
728          *  VARIABLE name
729          */
730         r = SDATA(root->this);
731 
732         /*
733          *  check for name conflicts
734          */
735         if (fnc_check(r) || com_check(r) || lst_find(CONSTANTS, r))
736         {
737           error( "VARIABLE not created [%s], identifier in use.\n ", r);
738         }
739 
740         if ((res = evaltree(LINK(root->this))) != NULL)
741         {
742           var_delete(r);
743           var = var_new(r,TYPE(res),1,1);
744 
745           d = MATR(res);
746           for(i = 0; i < NCOL(res)*NROW(res); i++)
747           {
748             *MATR(var) = *d++;
749             ptr = evalclause(LINK(root));
750           }
751 
752           var_delete_temp(res);
753         }
754         root = root->jmp;
755       break;
756       }
757     }
758     root = LINK(root);
759   }
760   return ptr;
761 }
762 
put_values(ptr,resname,par)763 VARIABLE *put_values(ptr, resname, par) VARIABLE *ptr, *par; char *resname;
764 /*======================================================================
765 ?  extract values from ptr, indexes in par, and put them to res
766 ^=====================================================================*/
767 {
768   static double defind = 0.0;
769 #pragma omp threadprivate (defind)
770 
771   double *ind1,
772          *ind2,
773          *dtmp;
774 
775   int i, j, k,
776       rows, cols,
777       size1, size2,
778       imax1, imax2,
779       ind;
780 
781   VARIABLE *res;
782 
783   res = (VARIABLE *)lst_find(VARIABLES, resname);
784 
785   if (NEXT(par) == NULL)
786   {
787     if (
788       res != NULL && NROW(par) == NROW(res) && NCOL(par) == NCOL(res)
789                   && !(NROW(res) == 1 && NCOL(res) == 1)
790       )
791     {
792       int logical = TRUE,
793           csize   = 0;
794 
795       dtmp = MATR(par);
796       for(i = 0; i < NROW(par)*NCOL(par); i++)
797         if (dtmp[i] != 0 && dtmp[i] != 1)
798         {
799           logical = FALSE;
800           break;
801         }
802 
803       if (logical)
804       {
805         imax1 = NROW(ptr) * NCOL(ptr);
806         dtmp  = MATR(ptr);
807         for(i = 0, k = 0; i < NROW(res); i++)
808           for(j = 0,csize=0; j < NCOL(res); j++)
809           {
810             while(M(par,i,j)==1 && j+csize<NCOL(res) && k+csize<imax1) csize++;
811             if (csize > 0)
812             {
813               memcpy(&M(res,i,j),&dtmp[k],csize*sizeof(double));
814               j += csize-1;
815               k += csize;
816               csize = 0;
817               if (k >= imax1) k = 0;
818             }
819           }
820          var_delete_temp(ptr);
821          return res;
822       }
823       else
824       {
825         ind1  = &defind;
826         size1 = 1;
827         ind2  = MATR(par);
828         size2 = NCOL(par);
829       }
830     }
831     else
832     {
833       ind1  = &defind;
834       size1 = 1;
835       ind2  = MATR(par);
836       size2 = NCOL(par);
837     }
838   }
839   else
840   {
841     ind1  = MATR(par);
842     size1 = NCOL(par);
843     ind2  = MATR(NEXT(par));
844     size2 = NCOL(NEXT(par));
845   }
846 
847   imax1 = (int)ind1[0];
848   for(i = 1; i < size1; i++)
849   {
850     imax1 = max(imax1, (int)ind1[i]);
851   }
852 
853   imax2 = (int)ind2[0];
854   for(i = 1; i < size2; i++)
855   {
856     imax2 = max(imax2, (int)ind2[i]);
857   }
858 
859   if (res == NULL)
860   {
861     res = var_new(resname, TYPE(ptr), imax1+1, imax2+1);
862   }
863   else if (NROW(res) <= imax1 || NCOL(res) <= imax2)
864   {
865     int ir = NROW(res), jc = NCOL(res);
866     MATRIX *t;
867 
868     imax1 = max(ir, imax1 + 1);
869     imax2 = max(jc, imax2 + 1);
870 
871     t = mat_new(TYPE(res), imax1, imax2);
872     dtmp = t->data;
873 
874     for(i = 0; i < ir; i++)
875     {
876       memcpy(&dtmp[i*imax2],&M(res,i,0),jc*sizeof(double));
877     }
878 
879     if (--REFCNT(res) == 0)
880     {
881       mat_free(res->this);
882     }
883     res->this = t;
884     REFCNT(res) = 1;
885   }
886   else if (REFCNT(res) > 1)
887   {
888     --REFCNT(res);
889     res->this = mat_copy(res->this);
890   }
891 
892   imax1 = NROW(ptr) * NCOL(ptr);
893   dtmp = MATR(ptr);
894   for(i = 0, k = 0; i < size1; i++)
895   {
896     ind = (int)ind1[i];
897     for(j = 0; j < size2; j++)
898     {
899       memcpy(&M(res,ind,(int)ind2[j]),&dtmp[k++],sizeof(double));
900       if (k >= imax1) k = 0;
901     }
902   }
903 
904   var_delete_temp(ptr);
905 
906   return res;
907 }
908 
put_result(ptr,resname,par,indexflag,printflag)909 VARIABLE *put_result(ptr, resname, par, indexflag, printflag)
910 /*======================================================================
911 ?  copy VARIABLE from one place to another
912 |  and conditionally print the result
913 ^=====================================================================*/
914    int indexflag, printflag;
915    VARIABLE *ptr, *par;
916    char *resname;
917 {
918    VARIABLE *res = NULL;
919 
920   var_delete( "ans" );
921   if (indexflag && par)
922   {
923     res = put_values(ptr, resname, par);
924   }
925   else
926   {
927     res = var_rename( ptr, resname );
928   }
929 
930   if ( res ) res->changed = 1;
931   if (printflag) var_print(res);
932 
933   return res;
934 }
935