1 #include <stddef.h>
2 #include <limits.h>
3
4 #include "ngspice/bool.h"
5 #include "ngspice/cpextern.h"
6 #include "ngspice/dvec.h"
7 #include "ngspice/fteext.h"
8 #include "ngspice/ngspice.h"
9 #include "ngspice/stringskip.h"
10
11 #include "com_display.h"
12 #include "com_let.h"
13 #include "completion.h"
14
15 /* Range of index values, such as 2:3 */
16 typedef struct index_range {
17 int low;
18 int high;
19 } index_range_t;
20
21 static void copy_vector_data(struct dvec *vec_dst,
22 const struct dvec *vec_src);
23 static void copy_vector_data_with_stride(struct dvec *vec_dst,
24 const struct dvec *vec_src,
25 int n_dst_index, const index_range_t *p_dst_index);
26 static int find_indices(char *s, const struct dvec *vec_dst,
27 index_range_t *p_index);
28 static int get_index_values(char *s, int n_elem_this_dim,
29 index_range_t *p_range);
30 int get_one_index_value(const char *s, int *p_index);
31
32 /* let <vec_name> = <expr>
33 * let <vec_name> = <vec_name_old> if variable 'plainlet' is set
34 * let <vec_name>[<bracket_expr>] = <expr>
35 * <bracket_expr> = <index_expr> <sep> <index_expr> <sep> ...
36 * <index_expr>
37 * <index_expr> = <expr> | <expr> : <expr>
38 * <sep> = "," | "] <ws> ["
39 * <expr> = standard ngspice expression
40 */
com_let(wordlist * wl)41 void com_let(wordlist *wl)
42 {
43 char *p, *s;
44 index_range_t p_dst_index[MAXDIMS];
45 int n_dst_index;
46 struct pnode *names = (struct pnode *) NULL;
47 struct dvec *vec_src = (struct dvec *) NULL;
48 char *rhs;
49
50 /* Start of index NULL is a flag for no index */
51 char *p_index_start = (char *) NULL;
52
53
54 /* let with no arguments is equivalent to display */
55 if (!wl) {
56 com_display(NULL);
57 return;
58 }
59
60 p = wl_flatten(wl); /* Everything after let -> string */
61
62 /* Separate vector name from RHS of assignment */
63 n_dst_index = 0;
64 if ((rhs = strchr(p, '=')) == (char *) NULL) {
65 fprintf(cp_err, "Error: bad let syntax\n");
66 txfree(p);
67 return;
68 }
69 *rhs++ = '\0';
70
71 /* Handle indexing. At start, p = LHS; rhs = RHS. If index is found
72 * p = leftmost part of orig p up to first '['. So p always
73 * becomes the vector name, possibly with some spaces at the end. */
74 if ((s = strchr(p, '[')) != NULL) {
75 /* This null makes the dest vector name a null-terminated string */
76 *s = '\0';
77 p_index_start = s + 1;
78 }
79
80 /* "Remove" any spaces at the end of the vector name at p by stepping
81 * from the end of the word to the first charatcter that is not
82 * whitespace and then overwriting the next character (which will be
83 * the original NULL if there was no whitespace) with a NULL. */
84 {
85 char *q;
86 for (q = p + strlen(p) - 1; *q <= ' ' && p <= q; q--) {
87 ;
88 }
89 *++q = '\0';
90 }
91
92 /* Sanity check */
93 if (eq(p, "all") || strchr(p, '@') || *p == '\0' || isdigit_c(*p)) {
94 fprintf(cp_err, "Error: bad variable name \"%s\"\n", p);
95 goto quit;
96 }
97
98 /* Locate the vector being assigned values. If NULL, the vector
99 * does not exist */
100 struct dvec *vec_dst = vec_get(p);
101 if (vec_dst != (struct dvec *) NULL) {
102 /* Fix-up dimension count and limit. Sometimes these are
103 * not set properly. If not set, give the vector 1 dimension and
104 * ensure the right length */
105 if (vec_dst->v_numdims < 1) {
106 vec_dst->v_numdims = 1;
107 }
108 if (vec_dst->v_numdims == 1) {
109 vec_dst->v_dims[0] = vec_dst->v_length;
110 }
111 }
112
113 /* If the vector was indexed, find the indices */
114 if (p_index_start != (char *) NULL) {
115 /* Test for an attempt to index an undefined vector */
116 if (vec_dst == (struct dvec *) NULL) {
117 fprintf(cp_err,
118 "When creating a new vector, it cannot be indexed.\n");
119 goto quit;
120 }
121
122 if (find_indices(p_index_start, vec_dst, p_dst_index) != 0) {
123 txfree(p);
124 return;
125 }
126 n_dst_index = vec_dst->v_numdims;
127 } /* end of case that an indexing bracket '[' was found */
128
129 /* Evaluate rhs */
130 /* Just copy a vector. rhs has to be a valid existing vector name
131 May be used to copy vectors with forbidden characters in their names
132 into a vector with a valid name.*/
133 if (cp_getvar("plainlet", CP_BOOL, NULL, 0)) {
134 vec_src = vec_get(rhs);
135 if (vec_src == (struct dvec *) NULL) {
136 fprintf(cp_err, "Error: Can't evaluate \"%s\"\n", rhs);
137 goto quit;
138 }
139 }
140 /* evaluate the rhs expression as usual, math characters may not be used in vec names,
141 the expression parser then will complain about a syntax error */
142 else {
143 if ((names = ft_getpnames_from_string(
144 rhs, TRUE)) == (struct pnode*)NULL) {
145 fprintf(cp_err, "Error: RHS \"%s\" invalid\n", rhs);
146 goto quit;
147 }
148 if ((vec_src = ft_evaluate(names)) == (struct dvec*)NULL) {
149 fprintf(cp_err, "Error: Can't evaluate \"%s\"\n", rhs);
150 goto quit;
151 }
152 }
153
154 if (vec_src->v_link2) {
155 fprintf(cp_err, "Warning: extra wildcard values ignored\n");
156 }
157
158 /* Fix-up dimension count and limit. Sometimes these are
159 * not set properly. If not set, give the vector 1 dimension and ensure
160 * the right length */
161 if (vec_src->v_numdims < 1) {
162 vec_src->v_numdims = 1;
163 }
164 if (vec_src->v_numdims == 1) {
165 vec_src->v_dims[0] = vec_src->v_length;
166 }
167
168 if (vec_dst == (struct dvec *) NULL) {
169 /* p is not an existing vector. So make a new one equal to vec_src
170 * in all ways, except enforce that it is a permanent vector. */
171 vec_dst = dvec_alloc(copy(p),
172 vec_src->v_type,
173 vec_src->v_flags | VF_PERMANENT,
174 vec_src->v_length, NULL);
175
176 copy_vector_data(vec_dst, vec_src);
177 vec_new(vec_dst); /* Add tp current plot */
178 cp_addkword(CT_VECTOR, vec_dst->v_name);
179 } /* end of case of new vector */
180 else {
181 /* Existing vector.*/
182
183 if (n_dst_index == 0) {
184 /* Not indexed, so make equal to source vector as if it
185 * was a new vector, except reuse the allocation if it
186 * is the same type (real/complex) and the allocation size
187 * is sufficient but not too large (>2X) . */
188 if (isreal(vec_dst) == isreal(vec_src) &&
189 vec_dst->v_alloc_length >= vec_src->v_length &&
190 vec_dst->v_alloc_length <= 2 * vec_src->v_length) {
191 vec_dst->v_length = vec_src->v_length;
192 copy_vector_data(vec_dst, vec_src);
193 }
194 else { /* Something not OK, so free and allocate again */
195 int n_elem_alloc = vec_src->v_alloc_length;
196 if (isreal(vec_dst)) {
197 tfree(vec_dst->v_realdata);
198 }
199 else { /* complex */
200 tfree(vec_dst->v_compdata);
201 }
202 if (isreal(vec_src)) {
203 vec_dst->v_realdata = TMALLOC(double, n_elem_alloc);
204 }
205 else { /* complex source */
206 vec_dst->v_compdata = TMALLOC(ngcomplex_t, n_elem_alloc);
207 }
208
209 /* Make the destination vector the right data type. A few
210 * extra () added to keep some compilers from warning. */
211 vec_dst->v_flags = (short int) (
212 ((int) vec_dst->v_flags & ~(VF_REAL | VF_COMPLEX)) |
213 ((int) vec_src->v_flags & (VF_REAL | VF_COMPLEX)));
214 vec_dst->v_alloc_length = vec_src->v_alloc_length;
215 vec_dst->v_length = vec_src->v_length;
216 copy_vector_data(vec_dst, vec_src);
217 }
218 }
219 /* Else indexed. In this case, the source data must fit the indexed
220 * range */
221 else {
222 {
223 int n_dst_elem = 1;
224 int i;
225 for (i = 0; i < n_dst_index; ++i) {
226 index_range_t *p_range_cur = p_dst_index + i;
227 n_dst_elem *= p_range_cur->high - p_range_cur->low + 1;
228 }
229
230 /* Check # elem required vs available */
231 if (n_dst_elem != vec_src->v_length) {
232 const int v_length = vec_src->v_length;
233 const bool f_1 = v_length == 1;
234 (void) fprintf(cp_err, "Data for an index vector must "
235 "fit exactly. The indexed range required %d "
236 "element%s to fill it, but there %s %d "
237 "element%s supplied.\n",
238 n_dst_elem, n_dst_elem == 1 ? "" : "s",
239 f_1 ? "was" : "were", v_length, f_1 ? "" : "s");
240 goto quit;
241 }
242 }
243
244 /* Real source data can be put into a complex destination,
245 * but the other way around is not possible */
246 if (isreal(vec_dst) && iscomplex(vec_src)) {
247 (void) fprintf(cp_err, "Complex data cannot be used "
248 "to fill an array of real data.\n");
249 goto quit;
250 }
251
252 /* OK to copy, so copy */
253 copy_vector_data_with_stride(vec_dst, vec_src,
254 n_dst_index, p_dst_index);
255 } /* end of indexed vector */
256 } /* end of existing vector */
257
258 vec_dst->v_minsignal = 0.0; /* How do these get reset ??? */
259 vec_dst->v_maxsignal = 0.0;
260 vec_dst->v_scale = vec_src->v_scale;
261
262 quit:
263 /* va: garbage collection for vec_src, if ft_evaluate() created a
264 * new vector while evaluating pnode `names' */
265 if (names != (struct pnode *) NULL) {
266 if (!names->pn_value && vec_src) {
267 vec_free(vec_src);
268 }
269 /* frees also vec_src, if pnode `names' is simple value */
270 free_pnode(names);
271 }
272 txfree(p);
273 } /* end of function com_let */
274
275
276
277 /* Process indexing portion of a let command. On entry, s is the address
278 * of the first byte after the first opening index bracket */
find_indices(char * s,const struct dvec * vec_dst,index_range_t * p_index)279 static int find_indices(char *s, const struct dvec *vec_dst,
280 index_range_t *p_index)
281 {
282 const int v_numdims_dst = vec_dst->v_numdims;
283 const int * const v_dims_dst = vec_dst->v_dims;
284 int dim_cur = 0; /* current dimension being set */
285
286
287 /* Can be either comma-separated or individual dimensions */
288 if (strchr(s, ',') != 0) { /* has commas */
289 char *p_end;
290 while ((p_end = strchr(s, ',')) != (char *) NULL) {
291 *p_end = '\0';
292 if (dim_cur == v_numdims_dst) {
293 (void) fprintf(cp_err, "Too many dimensions given.\n");
294 return -1;
295 }
296 if (get_index_values(s, v_dims_dst[dim_cur],
297 p_index + dim_cur) != 0) {
298 (void) fprintf(cp_err, "Dimension ranges "
299 "for dimension %d could not be found.\n",
300 dim_cur + 1);
301 return -1;
302 }
303 ++dim_cur;
304 s = p_end + 1; /* after (former) comma */
305 } /* end of loop over comma-separated indices */
306
307 /* Must be one more index ending with a bracket */
308 if ((p_end = strchr(s, ']')) == (char *) NULL) {
309 (void) fprintf(cp_err,
310 "Final dimension was not found.\n");
311 return -1;
312 }
313
314 *p_end = '\0';
315 if (dim_cur == v_numdims_dst) {
316 (void) fprintf(cp_err,
317 "Final dimension exceeded maximum number.\n");
318 return -1;
319 }
320 if (get_index_values(s, v_dims_dst[dim_cur],
321 p_index + dim_cur) != 0) {
322 (void) fprintf(cp_err, "Dimension ranges "
323 "for last dimension (%d) could not be found.\n",
324 dim_cur + 1);
325 return -1;
326 }
327 ++dim_cur;
328 s = p_end + 1;
329
330 /* Only white space is allowed after closing brace */
331 if (*(s = skip_ws(s)) != '\0') {
332 (void) fprintf(cp_err, "Invalid text was found "
333 "after dimension data for vector: \"%s\".\n",
334 s);
335 return -1;
336 }
337 } /* end of case x[ , , ] */
338 else { /* x[][][] */
339 char *p_end;
340 while ((p_end = strchr(s, ']')) != (char *) NULL) {
341 *p_end = '\0';
342 if (dim_cur == v_numdims_dst) {
343 (void) fprintf(cp_err, "Too many dimensions given. "
344 "Only %d are present.\n", v_numdims_dst);
345 return -1;
346 }
347 if (get_index_values(s, v_dims_dst[dim_cur],
348 p_index + dim_cur) != 0) {
349 (void) fprintf(cp_err, "Dimension ranges "
350 "for dimension %d could not be found.\n",
351 dim_cur + 1);
352 return -1;
353 }
354 ++dim_cur;
355 s = p_end + 1; /* after (former) ']' */
356 if (*(s = skip_ws(s)) == '\0') { /* reached end */
357 break;
358 }
359
360 /* Not end of expression, so must be '[' */
361 if (*s != '[') {
362 (void) fprintf(cp_err, "Dimension bracket '[' "
363 "for dimension %d could not be found.\n",
364 dim_cur + 1);
365 return -1;
366 }
367 s++; /* past '[' */
368 } /* end of loop over individual bracketed entries */
369
370 if (dim_cur == 0) {
371 /* Did not find a single ']' in the string */
372 (void) fprintf(cp_err, "The ']' for dimension 1 "
373 "could not be found.\n");
374 return -1;
375 }
376 } /* end of case x[][][][] */
377
378 /* Finalize dimensions. There must be the same number or one less than
379 * the number of dimensions of the vector. For the special case of one
380 * less, the final dimension is set to the full range of the last
381 * dimension. Note that checks for too many dimensions have already
382 * been performed while the index information was found. */
383 if (dim_cur != v_numdims_dst) { /* special case or error */
384 if (dim_cur == v_numdims_dst - 1) { /* special case */
385 index_range_t * const p_index_last = p_index + dim_cur;
386 p_index_last->low = 0;
387 p_index_last->high = vec_dst->v_dims[dim_cur] - 1;
388 }
389 else {
390 (void) fprintf(cp_err, "Error: Only %d dimensions "
391 "were supplied, but %d are needed. The last dimension "
392 "may be omitted, in which case it will default to the "
393 "full range of that dimension.\n",
394 dim_cur, v_numdims_dst);
395 }
396 }
397
398 return 0;
399 } /* end of function find_indices */
400
401
402
403 /* Convert expresion expr -> low and high ranges equal or
404 * expression expr1 : epr2 -> low = expr1 and high = expr2.
405 * Values are tested to ensure they are positive and that the low
406 * value does not exceed the high value.
407 *
408 * If expr1 is whitespace or empty, it defaults to 0. For high, the
409 * largest possible values is used.
410 */
get_index_values(char * s,int n_elem_this_dim,index_range_t * p_range)411 static int get_index_values(char *s, int n_elem_this_dim,
412 index_range_t *p_range)
413 {
414 char *p_colon;
415 if ((p_colon = strchr(s, ':')) == (char *) NULL) { /* One expression */
416 if (get_one_index_value(s, &p_range->low) != 0) {
417 (void) fprintf(cp_err, "Error getting index.\n");
418 return -1;
419 }
420 p_range->high = p_range->low;
421 }
422 else { /* l:h. If l defaults to 0 and h to dim size - 1 */
423 *p_colon = '\0';
424 {
425 const int rc = get_one_index_value(s, &p_range->low);
426 if (rc != 0) {
427 if (rc < 0) { /* error */
428 (void) fprintf(cp_err, "Error getting low range.\n");
429 return -1;
430 }
431 /* +1 -> Else use default */
432 p_range->low = 0;
433 }
434 }
435 s = p_colon + 1; /* past (former) colon */
436 {
437 const int rc = get_one_index_value(s, &p_range->high);
438 if (rc != 0) {
439 if (rc < 0) { /* error */
440 (void) fprintf(cp_err, "Error getting high range.\n");
441 return -1;
442 }
443 /* +1 -> Else use default */
444 p_range->high = n_elem_this_dim - 1;
445 }
446 }
447
448 /* Ensure ranges given were valid */
449 if (p_range->low > p_range->high) {
450 (void) fprintf(cp_err, "Error: low range (%d) is greater "
451 "than high range (%d).\n",
452 p_range->low, p_range->high);
453 return -1;
454 }
455 if (p_range->high >= n_elem_this_dim) {
456 (void) fprintf(cp_err, "Error: high range (%d) exceeds "
457 "the maximum value (%d).\n",
458 p_range->high, n_elem_this_dim - 1);
459 return -1;
460 }
461 }
462 return 0;
463 } /* end of function get_index_values */
464
465
466
467 /* Get an index value
468 *
469 * Return codes
470 * +1: String empty or all whitespace
471 * 0: Normal
472 * -1: Error
473 */
get_one_index_value(const char * s,int * p_index)474 int get_one_index_value(const char *s, int *p_index)
475 {
476 /* Test for a string of whitespace */
477 if (*(s = skip_ws(s)) == '\0') {
478 return +1;
479 }
480
481 /* Parse the expression */
482 struct pnode * const names = ft_getpnames_from_string(s, TRUE);
483 if (names == (struct pnode *) NULL) {
484 (void) fprintf(cp_err, "Unable to parse index expression.\n");
485 return -1;
486 }
487
488 /* Evaluate the parsing */
489 struct dvec * const t = ft_evaluate(names);
490 if (t == (struct dvec *) NULL) {
491 (void) fprintf(cp_err, "Unable to evaluate index expression.\n");
492 free_pnode_x(names);
493 return -1;
494 }
495
496 int xrc = 0;
497 if (t->v_link2 || t->v_length != 1 || !t->v_realdata) {
498 fprintf(cp_err, "Index expression is not a real scalar.\n");
499 xrc = -1;
500 }
501 else {
502 const int index = (int) floor(t->v_realdata[0] + 0.5);
503 if (index < 0) {
504 printf("Negative index (%d) is not allowed.\n", index);
505 xrc = -1;
506 }
507 else { /* index found ok */
508 *p_index = index;
509 }
510 }
511
512 /* Free resources */
513 if (names->pn_value != (struct dvec *) NULL) {
514 /* allocated value given to t */
515 vec_free_x(t);
516 }
517 free_pnode_x(names);
518
519 return xrc;
520 } /* end of function get_one_index_value */
521
522
523
524 /* Copy vector data and its metadata */
copy_vector_data(struct dvec * vec_dst,const struct dvec * vec_src)525 static void copy_vector_data(struct dvec *vec_dst,
526 const struct dvec *vec_src)
527 {
528 const size_t length = (size_t) vec_src->v_length;
529 int n_dim = vec_dst->v_numdims = vec_src->v_numdims;
530 (void) memcpy(vec_dst->v_dims, vec_src->v_dims,
531 (size_t) n_dim * sizeof(int));
532 if (isreal(vec_src)) {
533 (void) memcpy(vec_dst->v_realdata, vec_src->v_realdata,
534 length * sizeof(double));
535 }
536 else {
537 (void) memcpy(vec_dst->v_compdata, vec_src->v_compdata,
538 length * sizeof(ngcomplex_t));
539 }
540 } /* end of function copy_vector_data */
541
542
543
544 /* Copy vector data and its metadata using stride info */
copy_vector_data_with_stride(struct dvec * vec_dst,const struct dvec * vec_src,int n_dim,const index_range_t * p_range)545 static void copy_vector_data_with_stride(struct dvec *vec_dst,
546 const struct dvec *vec_src,
547 int n_dim, const index_range_t *p_range)
548 {
549 /* Offsets and related expressions at different levels of indexing
550 * given in elements
551 *
552 * Example
553 * Dimensions: 4
554 * Dimension extents: 10 X 8 X 100 X 5
555 * Selected ranges: 2:5 X 3:4 X 20:30 X 3:4
556 * Strides: 4000, 500, 5, 1
557 * Min offsets: 8000, 1500, 100, 3 -- offset to 1st
558 * element of range
559 * Cur cum offsets: 8000, 9500, 9600, 9603 (initial)
560 * Cur index: 2, 3, 20, X (initial)
561 *
562 * Note that the strides are built from the highest dimension,
563 * which always has stride 1, backwards.
564 */
565 int p_stride_level[MAXDIMS];
566 /* Stride changing index by 1 at each level */
567 int p_offset_level_min[MAXDIMS]; /* Offset to 1st elem at level */
568
569 /* Current cumulative offset at each level. A -1 index is created
570 * to handle the case of a single dimension more uniformly */
571 int p_offset_level_cum_full[MAXDIMS + 1];
572 int *p_offset_level_cum = p_offset_level_cum_full + 1;
573
574 int p_index_cur[MAXDIMS]; /* Current range value at each level */
575
576 {
577 const int index_max = n_dim - 1;
578 p_stride_level[index_max] = 1;
579 int *p_dim_ext = vec_dst->v_dims;
580 int i;
581 for (i = n_dim - 2; i >= 0; --i) {
582 const int i1 = i + 1;
583 p_stride_level[i] = p_stride_level[i1] * p_dim_ext[i1];
584 }
585 }
586
587 /* Initialize the minimum offsets, cumulative current offsets, and
588 * current index based on ranges and strides */
589 {
590 const int low_cur = p_index_cur[0] = p_range[0].low;
591 p_offset_level_cum[0] = p_offset_level_min[0] =
592 low_cur * p_stride_level[0];
593 }
594
595 {
596 int i;
597 for (i = 1; i < n_dim; ++i) {
598 const int low_cur = p_index_cur[i] = p_range[i].low;
599 p_offset_level_cum[i] = p_offset_level_cum[i - 1] +
600 (p_offset_level_min[i] = low_cur * p_stride_level[i]);
601 }
602 }
603
604 /* There are three cases to consider:
605 * 1) real dst <- real src
606 * 2) complex dst <- complex src
607 * 3) complex dst <- real src
608 *
609 * The first two can copy blocks at the highest dimesion and the can
610 * be combined by generalizing to the data size (sizeof(double) or
611 * sizeof(ngcomplex_t)) and offset of the data array. The third one
612 * must be assigned element by element with 0's given to the imaginary
613 * part of the data.
614 */
615
616 if (isreal(vec_src) && iscomplex(vec_dst)) {
617 /* complex dst <- real src */
618 int n_elem_topdim; /* # elements copied in top (stride 1) dimension */
619 ngcomplex_t *p_vec_data_dst = vec_dst->v_compdata;
620 /* Location of data in dvec struct */
621 double *p_vec_data_src = vec_src->v_realdata;
622 /* Location of data in dvec struct */
623
624 {
625 const int index_max = n_dim - 1;
626 const index_range_t * const p_range_max = p_range + index_max;
627 n_elem_topdim = p_range_max->high - p_range_max->low + 1;
628 }
629
630 /* Copy all data. Each loop iteration copies all of the elements
631 * at the highest dimension (which are contiguous). On entry to
632 * the loop, the arrays are initialized so that the first element
633 * can be copied, and they are updated in each iteration to
634 * process the next element. Note that if this function is called,
635 * there will always be at least one element to copy, so it
636 * is always safe to copy then check for the end of data. */
637 {
638 const int n_cpy = n_dim - 1; /* index where copying done */
639 const double *p_vec_data_src_end = p_vec_data_src +
640 vec_src->v_length; /* end of copying */
641 for ( ; ; ) {
642 /* Copy the data currently being located by the cumulative
643 * offset and the source location */
644 {
645 ngcomplex_t *p_dst_cur = p_vec_data_dst +
646 p_offset_level_cum[n_cpy];
647 ngcomplex_t *p_dst_end = p_dst_cur + n_elem_topdim;
648 for ( ; p_dst_cur < p_dst_end;
649 ++p_dst_cur, ++p_vec_data_src) {
650 p_dst_cur->cx_real = *p_vec_data_src;
651 p_dst_cur->cx_imag = 0.0;
652 }
653 }
654
655 /* Test for end of source data and exit if reached */
656 if (p_vec_data_src == p_vec_data_src_end) {
657 break; /* Copy is complete */
658 }
659
660 /* Move to the next destination location. Since the loop
661 * was not exited yet, it must exist */
662 {
663 int level_cur = n_cpy;
664
665 /* Move back to the first dimension that is not at its
666 * last element */
667 while (p_index_cur[level_cur] ==
668 p_range[level_cur].high) {
669 --level_cur;
670 }
671
672 /* Now at the first dimension level that is not full.
673 * Increment here and reset the highe ones to their
674 * minimum values to "count up." */
675 ++p_index_cur[level_cur];
676 p_offset_level_cum[level_cur] +=
677 p_stride_level[level_cur];
678 for (++level_cur; level_cur <= n_cpy; ++level_cur) {
679 p_index_cur[level_cur] = p_range[level_cur].low;
680 p_offset_level_cum[level_cur] =
681 p_offset_level_cum[level_cur - 1] +
682 p_offset_level_min[level_cur];
683 }
684 } /* end of block updating destination */
685 } /* end of loop copying from source to destination */
686 } /* end of block */
687 } /* end of case both real or complex */
688 else { /* Both real or complex (complex src and real dst not allowed) */
689 int n_byte_elem; /* Size of element */
690 int n_elem_topdim; /* # elements copied in top (stride 1) dimension */
691 int n_byte_topdim; /* contiguous bytes */
692 void *p_vec_data_dst; /* Location of data in dvec struct */
693 void *p_vec_data_src; /* Location of data in dvec struct */
694
695 {
696 const int index_max = n_dim - 1;
697 const index_range_t * const p_range_max = p_range + index_max;
698 n_elem_topdim = p_range_max->high - p_range_max->low + 1;
699 }
700
701 if (isreal(vec_src)) { /* Both real */
702 n_byte_elem = (int) sizeof(double);
703 n_byte_topdim = n_elem_topdim * (int) sizeof(double);
704 p_vec_data_dst = vec_dst->v_realdata;
705 p_vec_data_src = vec_src->v_realdata;
706 }
707 else {
708 n_byte_elem = (int) sizeof(ngcomplex_t);
709 n_byte_topdim = n_elem_topdim * (int) sizeof(ngcomplex_t);
710 p_vec_data_dst = vec_dst->v_compdata;
711 p_vec_data_src = vec_src->v_compdata;
712 }
713
714 /* Add the offset of the top dimension to all of the lower ones
715 * since it will always be added when copying */
716 {
717 int i;
718 const int n_max = n_dim - 1;
719 int offset_top = p_range[n_max].low;
720 p_offset_level_cum[-1] = offset_top;
721 for (i = 0; i < n_max; ++i) {
722 p_offset_level_cum[i] += offset_top;
723 }
724 }
725
726 /* Because the copies are being done in terms of bytes rather
727 * than complex data elements or real data elements, convert
728 * the strides and offsets from elements to bytes */
729 {
730 p_offset_level_cum[-1] *= n_byte_elem;
731 int i;
732 const int n_max = n_dim - 1;
733 for (i = 0; i < n_max; i++) {
734 p_stride_level[i] *= n_byte_elem;
735 p_offset_level_min[i] *= n_byte_elem;
736 p_offset_level_cum[i] *= n_byte_elem;
737 }
738 }
739
740 /* Copy all data. Each loop iteration copies all of the elements
741 * at the highest dimension (which are contiguous). On entry to
742 * the loop, the arrays are initialized so that the first element
743 * can be copied, and they are updated in each iteration to
744 * process the next element. Note that if this function is called,
745 * there will always be at least one element to copy, so it
746 * is always safe to copy then check for the end of data. */
747 {
748 const int n_cpy = n_dim - 2; /* index where copying done */
749 const void *p_vec_data_src_end = (char *) p_vec_data_src +
750 (size_t) (vec_src->v_length * n_byte_elem);
751 /* end of copying */
752 for ( ; ; ) {
753 /* Copy the data currently being located by the cumulative
754 * offset and the source location */
755 (void) memcpy(
756 (char *) p_vec_data_dst + p_offset_level_cum[n_cpy],
757 p_vec_data_src,
758 (size_t) n_byte_topdim);
759
760 /* Move to the next source data and exit the loop if
761 * the end is reached.
762 * NOTE: EXITING BEFORE UPDATING THE DESTINATION WILL
763 * PREVENT OVERRUNNING BUFFERS */
764 if ((p_vec_data_src = (char *) p_vec_data_src +
765 n_byte_topdim) == p_vec_data_src_end) {
766 break; /* Copy is complete */
767 }
768
769 /* Move to the next destination location. Since the loop
770 * was not exited yet, it must exist */
771 {
772 int level_cur = n_cpy;
773
774 /* Move back to the first dimension that is not at its
775 * last element */
776 while (p_index_cur[level_cur] ==
777 p_range[level_cur].high) {
778 --level_cur;
779 }
780
781 /* Now at the first dimension level that is not full.
782 * Increment here and reset the highe ones to their
783 * minimum values to "count up." */
784 ++p_index_cur[level_cur];
785 p_offset_level_cum[level_cur] +=
786 p_stride_level[level_cur];
787 for (++level_cur; level_cur <= n_cpy; ++level_cur) {
788 p_index_cur[level_cur] = p_range[level_cur].low;
789 p_offset_level_cum[level_cur] =
790 p_offset_level_cum[level_cur - 1] +
791 p_offset_level_min[level_cur];
792 }
793 } /* end of block updating destination */
794 } /* end of loop copying from source to destination */
795 } /* end of block */
796 } /* end of case both real or complex */
797 } /* end of function copy_vector_data_with_stride */
798
799
800
801