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