NDOSolver / FiOracle
Interfaces and Solvers for NonDifferentiable Optimization
Loading...
Searching...
No Matches
OPTvect.h
1/*--------------------------------------------------------------------------*/
2/*----------------------- File OPTvect.h -----------------------------------*/
3/*--------------------------------------------------------------------------*/
4/*-- --*/
5/*-- A bunch of little useful template inline functions. --*/
6/*-- --*/
7/*-- Scalar functions: --*/
8/*-- --*/
9/*-- ABS(), sgn(), min(), max(), Swap(), CeilDiv() --*/
10/*-- --*/
11/*-- Number-returning vector functions: --*/
12/*-- --*/
13/*-- Norm(), OneNorm(), INFNorm(), SumV() --*/
14/*-- MaxVecV(), MinVecV(), MaxVecI(), MinVecI() --*/
15/*-- ScalarProduct(), ScalarProduct[B[B]]() --*/
16/*-- --*/
17/*-- Vector operations: --*/
18/*-- --*/
19/*-- Vect[M]Assign[B[B]](), VectSum[B[B]](), VectSubtract[B[B]](), --*/
20/*-- Vect[I]Scale[B[B]](), VectAdd(), VectDiff(), VectMult(), --*/
21/*-- VectDivide(), VectXcg[B[B]]() --*/
22/*-- --*/
23/*-- Sparse/dense vector transformations: --*/
24/*-- --*/
25/*-- Sparsify(), SparsifyT(), SparsifyAT(), Densify(), Compact() --*/
26/*-- --*/
27/*-- Array manipulation: --*/
28/*-- --*/
29/*-- Merge(), ShiftVect(), RotateVect(), ShiftRVect(), RotateRVect() --*/
30/*-- --*/
31/*-- Searching: --*/
32/*-- --*/
33/*-- Match(), EqualVect() --*/
34/*-- BinSearch(), BinSearch1(), BinSearch2() --*/
35/*-- HeapIns(), HeapDel() --*/
36/*-- --*/
37/*-- Antonio Frangioni --*/
38/*-- Dipartimento di Informatica --*/
39/*-- Universita' di Pisa --*/
40/*-- --*/
41/*--------------------------------------------------------------------------*/
42/*--------------------------------------------------------------------------*/
43/*----------------------------- DEFINITIONS --------------------------------*/
44/*--------------------------------------------------------------------------*/
45/*--------------------------------------------------------------------------*/
46
47#ifndef __OPTvect
48 #define __OPTvect /* self-identification: #endif at the end of the file */
49
50/*--------------------------------------------------------------------------*/
51/*------------------------------ INCLUDES ----------------------------------*/
52/*--------------------------------------------------------------------------*/
53
54#include "OPTtypes.h"
55
56/*--------------------------------------------------------------------------*/
57/*--------------------------- NAMESPACE ------------------------------------*/
58/*--------------------------------------------------------------------------*/
59
60namespace OPTtypes_di_unipi_it
61{
62
63/*--------------------------------------------------------------------------*/
64/*----------------------------- FUNCTIONS ----------------------------------*/
65/*--------------------------------------------------------------------------*/
66/*-- --*/
67/*-- Scalar functions, like ABS(), min(), max(), sgn() ... --*/
68/*-- --*/
69/*--------------------------------------------------------------------------*/
70/*--------------------------------------------------------------------------*/
71
72template< class T >
73inline T ABS( const T x )
74{
75 return( x >= T( 0 ) ? x : -x );
76 }
77
78/*--------------------------------------------------------------------------*/
79
80template< class T >
81inline T sgn( const T x )
82{
83 return( x ? ( x > T( 0 ) ? T( 1 ) : T( -1 ) ) : T( 0 ) );
84 }
85
86/*--------------------------------------------------------------------------*/
87/* max() and min() are already defined on gcc 3.x */
88
89#ifndef __GNUC__
90
91template< class T >
92inline T min( const T x , const T y )
93{
94 return( x <= y ? x : y );
95 }
96
97template< class T >
98inline T max( const T x , const T y )
99{
100 return( x >= y ? x : y );
101 }
102
103#else
104#if __GNUC__ < 3
105
106template< class T >
107inline T min( const T x , const T y )
108{
109 return( x <= y ? x : y );
110 }
111
112template< class T >
113inline T max( const T x , const T y )
114{
115 return( x >= y ? x : y );
116 }
117
118#endif
119#endif
120
121/*--------------------------------------------------------------------------*/
122
123template< class T >
124inline void Swap( T &v1 , T &v2 )
125{
126 const T temp = v1;
127
128 v1 = v2;
129 v2 = temp;
130 }
131
132/*--------------------------------------------------------------------------*/
133
134template< class T1, class T2 >
135inline T1 CeilDiv( const T1 x , const T2 y )
136{
137 // returns the ceiling of (smallest integer number not smaller than) x / y,
138 // which have both to be integer types because the `%' operation is used
139
140 T1 temp = x / y;
141 if( x % y )
142 temp++;
143
144 return( temp );
145 }
146
147/*--------------------------------------------------------------------------*/
148/*--------------------------------------------------------------------------*/
149/*-- --*/
150/*-- Vector functions: norms, scalar products, assignments of vectors, --*/
151/*-- multiplication by a scalar and so on. The types involved in these --*/
152/*-- operations must support the elementary arithmetical operations '+', --*/
153/*-- '-', '*' and '/', as well as the assigment '='. --*/
154/*-- --*/
155/*-- Vectors are usually just a set of n (a parameter of the function) --*/
156/*-- consecutive elements of a certain type; however, in many cases it --*/
157/*-- is useful to "restrict" them to a subset of their entries. Hence, --*/
158/*-- (read-only) vectors of indices (cIndex_Set) play a special role; for --*/
159/*-- a vector g, with g{B} (where B is a vector of indices) we indicate --*/
160/*-- the "restricted" vector [ g[ B[ i ] ]. A typical reason for dealing --*/
161/*-- with "restricted" vectors is that "sparse" vectors, those having few --*/
162/*-- nonzeroes w.r.t. their length, are very common. --*/
163/*-- --*/
164/*-- Vector of indices are meant to be "infinity-terminated", i.e., an --*/
165/*-- Inf< Index >() must be found immediately after the last "valid" --*/
166/*-- Index, and sometimes they are required to be ordered, typically in --*/
167/*-- increasing sense. --*/
168/*-- --*/
169/*--------------------------------------------------------------------------*/
170/*--------------------------------------------------------------------------*/
171/*-- Number-returning functions (norms, scalar products ...) --*/
172/*--------------------------------------------------------------------------*/
173/*--------------------------------------------------------------------------*/
174
175template< class T >
176inline T Norm( const T *g , Index n )
177{
178 // returns Sum{ i = 0 .. n - 1 } g[ i ]^2 = g * g
179
180 T t = 0;
181 for( ; n-- ; ) {
182 const T tmp = *(g++);
183 t += tmp * tmp;
184 }
185
186 return( t );
187 }
188
189/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
190
191template< class T >
192inline T Norm( const T *g , cIndex_Set B )
193{
194 // Norm( g{B} )
195
196 T t = 0;
197 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; ) {
198 const T tmp = g[ h ];
199 t += tmp * tmp;
200 }
201
202 return( t );
203 }
204
205/*--------------------------------------------------------------------------*/
206
207template< class T >
208inline T OneNorm( const T *g , Index n )
209{
210 // returns Sum{ i = 0 .. n - 1 } ABS( g[ i ] )
211
212 T t = 0;
213 for( ; n-- ; )
214 t += ABS( *(g++) );
215
216 return( t );
217 }
218
219/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
220
221template< class T >
222inline T OneNorm( const T *g , cIndex_Set B )
223{
224 // OneNorm( g{B} )
225
226 T t = 0;
227 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
228 t += ABS( g[ h ] );
229
230 return( t );
231 }
232
233/*--------------------------------------------------------------------------*/
234
235template< class T >
236inline T INFNorm( const T *g , Index n )
237{
238 // returns Max{ i = 0 .. n - 1 } ABS( g[ i ] )
239
240 T t = 0;
241 for( ; n-- ; ) {
242 const T tmp = *(g++);
243 if( t < tmp )
244 t = tmp;
245 }
246
247 return( t );
248 }
249
250/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
251
252template< class T >
253inline T INFNorm( const T *g , cIndex_Set B )
254{
255 // INFNorm( g{B} )
256
257 T t = 0;
258 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; ) {
259 const T tmp = g[ h ];
260 if( t < tmp )
261 t = tmp;
262 }
263
264 return( t );
265 }
266
267/*--------------------------------------------------------------------------*/
268/*--------------------------------------------------------------------------*/
269
270template< class T >
271inline T SumV( const T *g , Index n )
272{
273 // returns Sum{ i = 0 .. n - 1 } g[ i ]
274
275 T t = 0;
276 for( ; n-- ; )
277 t += *(g++);
278
279 return( t );
280 }
281
282/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
283
284template< class T >
285inline T SumV( const T *g , cIndex_Set B )
286{
287 // SumV( g{B} )
288
289 T t = 0;
290 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
291 t += g[ h ];
292
293 return( t );
294 }
295
296/*--------------------------------------------------------------------------*/
297
298template< class T >
299inline T MaxVecV( const T *g , Index n )
300{
301 // returns Max{ i = 0 .. n - 1 } g[ i ]; n *must* be > 0
302
303 T max = *g;
304 for( ; --n ; )
305 if( *(++g) > max )
306 max = *g;
307
308 return( max );
309 }
310
311/*--------------------------------------------------------------------------*/
312
313template< class T >
314inline T MinVecV( const T *g , Index n )
315{
316 // returns Min{ i = 0 .. n - 1 } g[ i ]; n *must* be > 0
317
318 T min = *g;
319 for( ; --n ; )
320 if( *(++g) < min )
321 min = *g;
322
323 return( min );
324 }
325
326/*--------------------------------------------------------------------------*/
327
328template< class T >
329inline Index MaxVecI( const T *g , cIndex n )
330{
331 // returns the *index* of the maximum element of the n-vector v
332
333 T max = *g;
334 Index maxi = 0;
335 for( Index i = maxi ; ++i < n ; )
336 if( *(++g) > max )
337 {
338 maxi = i;
339 max = *g;
340 }
341
342 return( maxi );
343 }
344
345/*--------------------------------------------------------------------------*/
346
347template< class T >
348inline Index MinVecI( const T *g , cIndex n )
349{
350 // returns the *index* of the minimum element of the n-vector v
351
352 T min = *g;
353 Index mini = 0;
354 for( Index i = mini ; ++i < n ; )
355 if( *(++g) < max )
356 {
357 mini = i;
358 min = *g;
359 }
360
361 return( mini );
362 }
363
364/*--------------------------------------------------------------------------*/
365/*--------------------------------------------------------------------------*/
366
367template< class T1, class T2 >
368inline T1 ScalarProduct( const T1 *g1 , const T2 *g2 , Index n )
369{
370 // returns Sum{ i = 0 .. n - 1 } g1[ i ] * g2[ i ]
371
372 T1 t = 0;
373 for( ; n-- ; )
374 t += (*(g1++)) * (*(g2++));
375
376 return( t );
377 }
378
379/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
380
381template< class T1, class T2 >
382inline T1 ScalarProduct( const T1 *g1 , const T2 *g2 , cIndex_Set B )
383{
384 // ScalarProduct( g1 , g2{B} )
385
386 T1 t = 0;
387 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
388 t += (*(g1++)) * g2[ h ];
389
390 return( t );
391 }
392
393/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
394
395template< class T1, class T2 >
396inline T1 ScalarProduct( const T1 *g1 , cIndex_Set B1 ,
397 const T2 *g2 , cIndex_Set B2 )
398{
399 // ScalarProduct( g1{B} , g2{B} ), B = intersection of B1 and B2
400
401 T1 t = 0;
402 Index h = *B1;
403 Index k = *B2;
404 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
405 for(;;)
406 if( k < h )
407 {
408 if( ( k = *(++B2) ) == Inf< Index >() )
409 break;
410 g2++;
411 }
412 else
413 if( h < k )
414 {
415 if( ( h = *(++B1) ) == Inf< Index >() )
416 break;
417 g1++;
418 }
419 else
420 {
421 t += (*(g1++)) * (*(g2++));
422 if( ( k = *(++B2) ) == Inf< Index >() )
423 break;
424 if( ( h = *(++B1) ) == Inf< Index >() )
425 break;
426 }
427 }
428
429 return( t );
430
431 } // end( ScalarProduct( g1 , B1 , g2 , B2 )
432
433/*--------------------------------------------------------------------------*/
434/*--------------------------------------------------------------------------*/
435
436template< class T1, class T2 >
437inline T1 ScalarProductB( const T1 *g1 , const T2 *g2 , cIndex_Set B )
438{
439 // ScalarProduct( g1{B} , g2 )
440
441 T1 t = 0;
442 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
443 t += g1[ h ] * (*(g2++));
444
445 return( t );
446 }
447
448/*--------------------------------------------------------------------------*/
449/*--------------------------------------------------------------------------*/
450
451template< class T1, class T2 >
452inline T1 ScalarProductBB( const T1 *g1 , const T2 *g2 , cIndex_Set B )
453{
454 // ScalarProduct( g1{B} , g2{B} )
455
456 T1 t = 0;
457 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
458 t += g1[ h ] * g2[ h ];
459
460 return( t );
461 }
462
463/*--------------------------------------------------------------------------*/
464/*--------------------------------------------------------------------------*/
465/*-- Vector operations: sum/difference of vectors, multiplication by a --*/
466/*-- scalar, assignment, setting, scaling ... --*/
467/*--------------------------------------------------------------------------*/
468/*--------------------------------------------------------------------------*/
469
470template< class T >
471inline void VectAssign( T *const g , const T x , cIndex n )
472{
473 // g[ i ] = x for each i = 0 .. n - 1
474
475 for( T *tg = g + n ; tg > g ; )
476 *(--tg) = x;
477 }
478
479/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
480
481template< class T >
482inline void VectAssign( T *const g , const T x , cIndex_Set B )
483{
484 // g{B} = x, all other entries of g unchanged
485
486 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
487 g[ h ] = x;
488 }
489
490/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
491
492template< class T1, class T2 >
493inline void VectAssign( T1 *g1 , const T2 *g2 , Index n )
494{
495 // g1 := g2
496
497 for( ; n-- ; )
498 *(g1++) = *(g2++);
499 }
500
501/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
502
503inline Index_Set VectAssign( Index_Set g1 , cIndex_Set g2 )
504{
505 // special version of the above for g2 an Inf< Index >()-terminated vector of
506 // indices do not write the terminating Inf< Index >(), but returns the
507 // pointer to the position where it should be written
508
509 for( Index h ; ( h = *(g2++) ) < Inf< Index >() ; )
510 *(g1++) = h;
511
512 return( g1 );
513 }
514
515/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
516
517template< class T1, class T2 >
518inline void VectAssign( T1 *g1 , const T2 *g2 , cIndex_Set B )
519{
520 // g1 = g2{B}
521
522 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
523 *(g1++) = g2[ h ];
524 }
525
526/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
527
528template< class T1, class T2 >
529inline void VectAssign( T1 *g1 , cIndex_Set B1 , const T2 *g2 , cIndex_Set B2 )
530{
531 // g1{B} = g2{B}, B = intersection of B1 and B2, all other unchanged
532
533 Index h = *B1;
534 Index k = *B2;
535 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
536 for(;;)
537 if( k < h )
538 {
539 if( ( k = *(++B2) ) == Inf< Index >() )
540 break;
541 g2++;
542 }
543 else
544 if( h < k )
545 {
546 if( ( h = *(++B1) ) == Inf< Index >() )
547 break;
548 g1++;
549 }
550 else
551 {
552 (*(g1++)) = (*(g2++));
553 if( ( k = *(++B2) ) == Inf< Index >() )
554 break;
555 if( ( h = *(++B1) ) == Inf< Index >() )
556 break;
557 }
558 }
559 } // end( VectAssign( g1 , B1 , g2 , B2 )
560
561/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
562
563template< class T1, class T2 >
564inline void VectAssign( T1 *g1 , cIndex_Set B1 , const T2 *g2 ,
565 cIndex_Set B2 , const T1 gg )
566{
567 // g1{B} = g2{B}, B = intersection of B1 and B2, g1{B1 / B} = gg
568
569 Index k = *B2;
570 for( Index h ; ( h = *(B1++) ) < Inf< Index >() ; ) {
571 while( k < h ) {
572 k = *(++B2);
573 g2++;
574 }
575
576 if( k == Inf< Index >() )
577 break;
578
579 if( h == k )
580 *(g1++) = *(g2++);
581 else
582 *(g1++) = gg;
583 }
584 } // end( VectAssign( g1 , B1 , g2 , B2 , gg )
585
586/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
587
588template< class T1 , class T2 , class T3 >
589inline void VectAssign( T1 *g1 , const T2 *g2 , const T3 x , Index n )
590{
591 // g1 := x * g2
592
593 for( ; n-- ; )
594 *(g1++) = x * (*(g2++));
595 }
596
597/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
598
599template< class T1 , class T2 , class T3 >
600inline void VectAssign( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
601{
602 // g1 = x * g2{B}
603
604 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
605 *(g1++) = x * g2[ h ];
606 }
607
608/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
609
610template< class T1 , class T2 , class T3 >
611inline void VectAssign( T1 *g1 , cIndex_Set B1 , const T3 x , const T2 *g2 ,
612 cIndex_Set B2 )
613{
614 // g1{B} = x * g2{B}, B = intersection of B1 and B2, all other unchanged
615
616 Index h = *B1;
617 Index k = *B2;
618 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
619 for(;;)
620 if( k < h ) {
621 if( ( k = *(++B2) ) == Inf< Index >() )
622 break;
623 g2++;
624 }
625 else
626 if( h < k ) {
627 if( ( h = *(++B1) ) == Inf< Index >() )
628 break;
629 g1++;
630 }
631 else {
632 (*(g1++)) = x * (*(g2++));
633 if( ( k = *(++B2) ) == Inf< Index >() )
634 break;
635 if( ( h = *(++B1) ) == Inf< Index >() )
636 break;
637 }
638 }
639 } // end( VectAssign( g1 , B1 , x , g2 , B2 )
640
641/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
642
643template< class T1 , class T2 , class T3 >
644inline void VectAssign( T1 *g1 , cIndex_Set B1 , const T3 x , const T2 *g2 ,
645 cIndex_Set B2 , const T1 gg )
646{
647 // g1{B} = x * g2{B}, B = intersection of B1 and B2, g1{B1 / B} = gg
648
649 Index k = *B2;
650 for( Index h ; ( h = *(B1++) ) < Inf< Index >() ; ) {
651 while( k < h ) {
652 k = *(++B2);
653 g2++;
654 }
655
656 if( k == Inf< Index >() )
657 break;
658
659 if( h == k )
660 *(g1++) = x * (*(g2++));
661 else
662 *(g1++) = gg;
663 }
664 } // end( VectAssign( g1 , B1 , x , g2 , B2 , gg )
665
666/*--------------------------------------------------------------------------*/
667
668template< class T >
669inline void VectAssignB( T *const g , const T x , cIndex_Set B )
670{
671 // g{B} = x, all other entries of g unchanged
672
673 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
674 g[ h ] = x;
675 }
676
677/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
678
679template< class T >
680inline void VectAssignB( T *g , const T x , cIndex_Set B , cIndex n ,
681 const T gg = 0 )
682{
683 // g{B} = x, all other entries (0 .. n - 1) of g1 = gg
684
685 Index h = *(B++);
686 for( Index i = 0 ; i < n ; )
687 if( h == i++ )
688 {
689 *(g++) = x;
690 h = *(B++);
691 }
692 else
693 *(g++) = gg;
694 }
695
696/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
697
698template< class T1, class T2 >
699inline void VectAssignB( T1 *g1 , const T2 *g2 , cIndex_Set B )
700{
701 // g1{B} = g2, all other entries of g1 unchanged
702
703 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
704 g1[ h ] = *(g2++);
705 }
706
707/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
708
709template< class T1, class T2 >
710inline void VectAssignB( T1 *g1 , const T2 *g2 , cIndex_Set B , cIndex n ,
711 const T1 gg = 0 )
712{
713 // g1{B} = g2, all other entries (0 .. n - 1) of g1 = gg
714
715 Index h = *(B++);
716 for( Index i = 0 ; i < n ; )
717 if( h == i++ ) {
718 *(g1++) = *(g2++);
719 h = *(B++);
720 }
721 else
722 *(g1++) = gg;
723 }
724
725/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
726
727template< class T1, class T2 >
728inline void VectAssignB( T1 *g1 , const T2 *g2 ,
729 const T1 x , cIndex_Set B )
730{
731 // g1{B} = g2, all other entries of g1 = unchanged
732
733 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
734 g1[ h ] = x * (*(g2++));
735 }
736
737/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
738
739template< class T1, class T2 >
740inline void VectAssignB( T1 *g1 , const T2 *g2 , const T1 x , cIndex_Set B ,
741 cIndex n , const T1 gg = 0 )
742{
743 // g1{B} = x * g2, all other entries (0 .. n - 1) of g1 = gg
744
745 Index h = *(B++);
746 for( Index i = 0 ; i < n ; )
747 if( h == i++ )
748 {
749 *(g1++) = x * (*(g2++));
750 h = *(B++);
751 }
752 else
753 *(g1++) = gg;
754 }
755
756/*--------------------------------------------------------------------------*/
757
758template< class T1, class T2 >
759inline void VectAssignBB( T1 *g1 , const T2 *g2 , cIndex_Set B )
760{
761 // g1{B} = g2{B}, all other entries unchanged
762
763 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
764 g1[ h ] = g2[ h ];
765 }
766
767/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
768
769template< class T1, class T2 >
770inline void VectAssignBB( T1 *g1 , const T2 *g2 , const T1 x , cIndex_Set B )
771{
772 // g1{B} = x * g2{B}, all other entries unchanged
773
774 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
775 g1[ h ] = x * g2[ h ];
776 }
777
778/*--------------------------------------------------------------------------*/
779/*--------------------------------------------------------------------------*/
780
781template< class T1, class T2 >
782inline void VectMAssign( T1 *g1 , const T2 *g2 , Index n )
783{
784 // g1 := - g2
785
786 for( ; n-- ; )
787 *(g1++) = - *(g2++);
788 }
789
790/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
791
792template< class T1, class T2 >
793inline void VectMAssign( T1 *g1 , const T2 *g2 , cIndex_Set B )
794{
795 // g1 = - g2{B}
796
797 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
798 *(g1++) = - g2[ h ];
799 }
800
801/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
802
803template< class T1, class T2 >
804inline void VectMAssign( T1 *g1 , cIndex_Set B1 , T2 *g2 , cIndex_Set B2 )
805{
806 // g1{B} = - g2{B}, B = intersection of B1 and B2, all other unchanged
807
808 Index h = *B1;
809 Index k = *B2;
810 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
811 for(;;)
812 if( k < h ) {
813 if( ( k = *(++B2) ) == Inf< Index >() )
814 break;
815 g2++;
816 }
817 else
818 if( h < k ) {
819 if( ( h = *(++B1) ) == Inf< Index >() )
820 break;
821 g1++;
822 }
823 else {
824 (*(g1++)) = - (*(g2++));
825 if( ( k = *(++B2) ) == Inf< Index >() )
826 break;
827 if( ( h = *(++B1) ) == Inf< Index >() )
828 break;
829 }
830 }
831 } // end( VectMAssign( g1 , B1 , g2 , B2 )
832
833/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
834
835template< class T1, class T2 >
836inline void VectMAssign( T1 *g1 , cIndex_Set B1 , const T2 *g2 ,
837 cIndex_Set B2 , const T1 gg )
838{
839 // g1{B} = - g2{B}, B = intersection of B1 and B2, g1{B1 / B} = gg
840
841 Index k = *B2;
842 for( Index h ; ( h = *(B1++) ) < Inf< Index >() ; ) {
843 while( k < h ) {
844 k = *(++B2);
845 g2++;
846 }
847
848 if( k == Inf< Index >() )
849 break;
850
851 if( h == k )
852 *(g1++) = - *(g2++);
853 else
854 *(g1++) = gg;
855 }
856 } // end( VectMAssign( g1 , B1 , g2 , B2 , gg )
857
858/*--------------------------------------------------------------------------*/
859
860template< class T1, class T2 >
861inline void VectMAssignB( T1 *g1 , const T2 *g2 , cIndex_Set B )
862{
863 // g1{B} = - g2, all other entries of g1 = unchanged
864
865 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
866 *(g1++) = - g2[ h ];
867 }
868
869/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
870
871template< class T1, class T2 >
872inline void VectMAssignB( T1 *g1 , const T2 *g2 , cIndex_Set B , cIndex n ,
873 const T1 gg = 0 )
874{
875 // g1{B} = - g2, all other entries (0 .. n - 1) of g1 = gg
876
877 Index h = *(B++);
878 for( Index i = 0 ; i < n ; )
879 if( h == i++ ) {
880 *(g1++) = - *(g2++);
881 h = *(B++);
882 }
883 else
884 *(g1++) = gg;
885 }
886
887/*--------------------------------------------------------------------------*/
888
889template< class T1, class T2 >
890inline void VectMAssignBB( T1 *g1 , const T2 *g2 , cIndex_Set B )
891{
892 // g1{B} = - g2{B}, all other entries unchanged
893
894 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
895 g1[ h ] = - g2[ h ];
896 }
897
898/*--------------------------------------------------------------------------*/
899/*--------------------------------------------------------------------------*/
900
901template< class T >
902inline void VectSum( T *const g , const T x , cIndex n )
903{
904 // g[ i ] += x for each i = 0 .. n - 1
905
906 for( T *tg = g + n ; tg > g ; )
907 *(--tg) += x;
908 }
909
910/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
911
912inline void VectSum( Index_Set g , cIndex x )
913{
914 // special version of the above for g an Inf< Index >()-terminated vector of
915 // indices
916
917 while( *g < Inf< Index >() )
918 *(g++) += x;
919 }
920
921/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
922
923template< class T >
924inline void VectSum( T *const g , const T x , cIndex_Set B )
925{
926 // g{B} += x, all other entries of g unchanged
927
928 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
929 g[ h ] += x;
930 }
931
932/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
933
934template< class T1, class T2 >
935inline void VectSum( T1 *g1 , const T2 *g2 , Index n )
936{
937 // g1 += g2
938
939 for( ; n-- ; )
940 *(g1++) += *(g2++);
941 }
942
943/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
944
945template< class T1, class T2 >
946inline void VectSum( T1 *g1 , const T2 *g2 , cIndex_Set B )
947{
948 // g1 += g2{B}
949
950 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
951 *(g1++) += g2[ h ];
952 }
953
954/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
955
956template< class T1, class T2 >
957inline void VectSum( T1 *g1 , cIndex_Set B1 , const T2 *g2 , cIndex_Set B2 )
958{
959 // g1{B} += g2{B}, B = intersection of B1 and B2, all other unchanged
960
961 Index h = *B1;
962 Index k = *B2;
963 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
964 for(;;)
965 if( k < h ) {
966 if( ( k = *(++B2) ) == Inf< Index >() )
967 break;
968 g2++;
969 }
970 else
971 if( h < k ) {
972 if( ( h = *(++B1) ) == Inf< Index >() )
973 break;
974 g1++;
975 }
976 else {
977 (*(g1++)) += (*(g2++));
978 if( ( k = *(++B2) ) == Inf< Index >() )
979 break;
980 if( ( h = *(++B1) ) == Inf< Index >() )
981 break;
982 }
983 }
984 } // end( VectSum( g1 , B1 , g2 , B2 )
985
986/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
987
988template< class T1 , class T2 , class T3 >
989inline void VectSum( T1 *g1 , const T2 *g2 , const T3 x , Index n )
990{
991 // g1 += x * g2
992
993 for( ; n-- ; )
994 *(g1++) += x * (*(g2++));
995 }
996
997/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
998
999template< class T1 , class T2 , class T3 >
1000inline void VectSum( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1001{
1002 // g1 += x * g2{B}
1003
1004 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1005 *(g1++) += g2[ h ];
1006 }
1007
1008/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1009
1010template< class T1 , class T2 , class T3 >
1011inline void VectSum( T1 *g1 , cIndex_Set B1 , const T3 x , const T2 *g2 ,
1012 cIndex_Set B2 )
1013{
1014 // g1{B} += x * g2{B}, B = intersection of B1 and B2, all other unchanged
1015
1016 Index h = *B1;
1017 Index k = *B2;
1018 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
1019 for(;;)
1020 if( k < h ) {
1021 if( ( k = *(++B2) ) == Inf< Index >() )
1022 break;
1023 g2++;
1024 }
1025 else
1026 if( h < k ) {
1027 if( ( h = *(++B1) ) == Inf< Index >() )
1028 break;
1029 g1++;
1030 }
1031 else {
1032 (*(g1++)) += x * (*(g2++));
1033 if( ( k = *(++B2) ) == Inf< Index >() )
1034 break;
1035 if( ( h = *(++B1) ) == Inf< Index >() )
1036 break;
1037 }
1038 }
1039 } // end( VectSum( g1 , B1 , x , g2 , B2 )
1040
1041/*--------------------------------------------------------------------------*/
1042
1043template< class T >
1044inline void VectSumB( T *g , const T x , cIndex_Set B , cIndex n , const T gg )
1045{
1046 // g{B} += x, all other entries (0 .. n - 1) of g1 += gg
1047
1048 Index h = *(B++);
1049 for( Index i = 0 ; i < n ; )
1050 if( h == i++ ) {
1051 *(g++) += x;
1052 h = *(B++);
1053 }
1054 else
1055 *(g++) += gg;
1056 }
1057
1058/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1059
1060template< class T1, class T2 >
1061inline void VectSumB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1062{
1063 // g1{B} += g2, all other entries of g1 unchanged
1064
1065 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1066 g1[ h ] += *(g2++);
1067 }
1068
1069/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1070
1071template< class T1, class T2 >
1072inline void VectSumB( T1 *g1 , const T2 *g2 , cIndex_Set B , cIndex n ,
1073 const T1 gg )
1074{
1075 // g1{B} += g2, all other entries (0 .. n - 1) of g1 += gg
1076
1077 Index h = *(B++);
1078 for( Index i = 0 ; i < n ; )
1079 if( h == i++ ) {
1080 *(g1++) += *(g2++);
1081 h = *(B++);
1082 }
1083 else
1084 *(g1++) += gg;
1085 }
1086
1087/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1088
1089template< class T1 , class T2 , class T3 >
1090inline void VectSumB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1091{
1092 // g1{B} += x * g2, all other entries of g1 unchanged
1093
1094 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1095 g1[ h ] += x * (*(g2++));
1096 }
1097
1098/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1099
1100template< class T1 , class T2 , class T3 >
1101inline void VectSumB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B ,
1102 cIndex n , const T1 gg )
1103{
1104 // g1{B} += x * g2, all other entries (0 .. n - 1) of g1 += gg
1105
1106 Index h = *(B++);
1107 for( Index i = 0 ; i < n ; )
1108 if( h == i++ ) {
1109 *(g1++) += x * (*(g2++));
1110 h = *(B++);
1111 }
1112 else
1113 *(g1++) += gg;
1114 }
1115
1116/*--------------------------------------------------------------------------*/
1117
1118template< class T1, class T2 >
1119inline void VectSumBB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1120{
1121 // g1{B} += g2{B}, all other entries unchanged
1122
1123 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1124 g1[ h ] += g2[ h ];
1125 }
1126
1127/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1128
1129template< class T1 , class T2 , class T3 >
1130inline void VectSumBB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1131{
1132 // g1{B} += x * g2{B}, all other entries unchanged
1133
1134 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1135 g1[ h ] += x * g2[ h ];
1136 }
1137
1138/*--------------------------------------------------------------------------*/
1139/*--------------------------------------------------------------------------*/
1140
1141template< class T >
1142inline void VectSubtract( T *const g , const T x , cIndex n )
1143{
1144 // g[ i ] -= x for each i = 0 .. n - 1; useful for *unsigned* data types
1145 // where VectSum( g , -x , n ) would not work
1146
1147 for( T *tg = g + n ; tg > g ; )
1148 *(--tg) -= x;
1149 }
1150
1151/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1152
1153inline void VectSubtract( Index_Set g , cIndex x )
1154{
1155 // special version of the above for g an Inf< Index >()-terminated vector of
1156 // indices (note that indices are generally unsigned)
1157
1158 while( *g < Inf< Index >() )
1159 *(g++) -= x;
1160 }
1161
1162/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1163
1164template< class T >
1165inline void VectSubtract( T *const g , const T x , cIndex_Set B )
1166{
1167 // g{B} -= x, all other entries of g unchanged
1168
1169 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1170 g[ h ] -= x;
1171 }
1172
1173/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1174
1175template< class T1, class T2 >
1176inline void VectSubtract( T1 *g1 , const T2 *g2 , Index n )
1177{
1178 // g1 -= g2 (element-wise)
1179
1180 for( ; n-- ; )
1181 *(g1++) -= *(g2++);
1182 }
1183
1184/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1185
1186template< class T1, class T2 >
1187inline void VectSubtract( T1 *g1 , const T2 *g2 , cIndex_Set B )
1188{
1189 // g1 -= g2{B} (element-wise)
1190
1191 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1192 *(g1++) -= g2[ h ];
1193 }
1194
1195/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1196
1197template< class T1, class T2 >
1198inline void VectSubtract( T1 *g1 , cIndex_Set B1 , const T2 *g2 ,
1199 cIndex_Set B2 )
1200{
1201 // g1{B} -= g2{B}, B = intersection of B1 and B2, all other unchanged
1202
1203 Index h = *B1;
1204 Index k = *B2;
1205 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
1206 for(;;)
1207 if( k < h ) {
1208 if( ( k = *(++B2) ) == Inf< Index >() )
1209 break;
1210 g2++;
1211 }
1212 else
1213 if( h < k ) {
1214 if( ( h = *(++B1) ) == Inf< Index >() )
1215 break;
1216 g1++;
1217 }
1218 else {
1219 (*(g1++)) -= (*(g2++));
1220 if( ( k = *(++B2) ) == Inf< Index >() )
1221 break;
1222 if( ( h = *(++B1) ) == Inf< Index >() )
1223 break;
1224 }
1225 }
1226 } // end( VectSubtract( g1 , B1 , g2 , B2 )
1227
1228/*--------------------------------------------------------------------------*/
1229
1230template< class T >
1231inline void VectSubtractB( T *g , const T x , cIndex_Set B , cIndex n ,
1232 const T gg )
1233{
1234 // g{B} -= x, all other entries (0 .. n - 1) of g1 -= gg
1235
1236 Index h = *(B++);
1237 for( Index i = 0 ; i < n ; )
1238 if( h == i++ ) {
1239 *(g++) -= x;
1240 h = *(B++);
1241 }
1242 else
1243 *(g++) -= gg;
1244 }
1245
1246/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1247
1248template< class T1, class T2 >
1249inline void VectSubtractB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1250{
1251 // g1{B} -= g2, all other entries of g1 unchanged
1252
1253 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1254 g1[ h ] -= *(g2++);
1255 }
1256
1257/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1258
1259template< class T1, class T2 >
1260inline void VectSubtractB( T1 *g1 , const T2 *g2 , cIndex_Set B , cIndex n ,
1261 const T1 gg )
1262{
1263 // g1{B} -= g2, all other entries (0 .. n - 1) of g1 -= gg
1264
1265 Index h = *(B++);
1266 for( Index i = 0 ; i < n ; )
1267 if( h == i++ ) {
1268 *(g1++) -= *(g2++);
1269 h = *(B++);
1270 }
1271 else
1272 *(g1++) -= gg;
1273 }
1274
1275/*--------------------------------------------------------------------------*/
1276
1277template< class T1, class T2 >
1278inline void VectSubtractBB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1279{
1280 // g1{B} -= g2{B}, all other entries unchanged
1281
1282 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1283 g1[ h ] -= g2[ h ];
1284 }
1285
1286/*--------------------------------------------------------------------------*/
1287/*--------------------------------------------------------------------------*/
1288
1289template< class T >
1290inline void VectScale( T *g , const T x , Index n )
1291{
1292 // g *= x, x a scalar
1293
1294 for( ; n-- ; )
1295 *(g++) *= x;
1296 }
1297
1298/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1299
1300template< class T >
1301inline void VectScale( T *g , const T x , cIndex_Set B )
1302{
1303 // g{B} *= x, x a scalar
1304
1305 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1306 g[ h ] *= x;
1307 }
1308
1309/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1310
1311template< class T1, class T2 >
1312inline void VectScale( T1 *g1 , const T2 *g2 , Index n )
1313{
1314 // g1 *= g2
1315
1316 for( ; n-- ; )
1317 *(g1++) *= *(g2++);
1318 }
1319
1320/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1321
1322template< class T1, class T2 >
1323inline void VectScale( T1 *g1 , const T2 *g2 , cIndex_Set B )
1324{
1325 // g1 *= g2{B}
1326
1327 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1328 *(g1++) *= g2[ h ];
1329 }
1330
1331/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1332
1333template< class T1, class T2 >
1334inline void VectScale( T1 *g1 , cIndex_Set B1 , const T2 *g2 , cIndex_Set B2 )
1335{
1336 // g1{B} *= g2{B}, B = intersection of B1 and B2, all other unchanged
1337
1338 Index h = *B1;
1339 Index k = *B2;
1340 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
1341 for(;;)
1342 if( k < h ) {
1343 if( ( k = *(++B2) ) == Inf< Index >() )
1344 break;
1345 g2++;
1346 }
1347 else
1348 if( h < k ) {
1349 if( ( h = *(++B1) ) == Inf< Index >() )
1350 break;
1351 g1++;
1352 }
1353 else {
1354 (*(g1++)) *= (*(g2++));
1355 if( ( k = *(++B2) ) == Inf< Index >() )
1356 break;
1357 if( ( h = *(++B1) ) == Inf< Index >() )
1358 break;
1359 }
1360 }
1361 } // end( VectScale( g1 , B1 , g2 , B2 )
1362
1363/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1364
1365template< class T1 , class T2 , class T3 >
1366inline void VectScale( T1 *g1 , const T2 *g2 , const T3 x , Index n )
1367{
1368 // g1 *= x * g2
1369
1370 for( ; n-- ; )
1371 *(g1++) *= x * (*(g2++));
1372 }
1373
1374/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1375
1376template< class T1 , class T2 , class T3 >
1377inline void VectScale( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1378{
1379 // g1 *= x * g2{B}
1380
1381 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1382 *(g1++) *= g2[ h ];
1383 }
1384
1385/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1386
1387template< class T1 , class T2 , class T3 >
1388inline void VectScale( T1 *g1 , cIndex_Set B1 , const T3 x , const T2 *g2 ,
1389 cIndex_Set B2 )
1390{
1391 // g1{B} *= x * g2{B}, B = intersection of B1 and B2, all other unchanged
1392
1393 Index h = *B1;
1394 Index k = *B2;
1395 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
1396 for(;;)
1397 if( k < h ) {
1398 if( ( k = *(++B2) ) == Inf< Index >() )
1399 break;
1400 g2++;
1401 }
1402 else
1403 if( h < k ) {
1404 if( ( h = *(++B1) ) == Inf< Index >() )
1405 break;
1406 g1++;
1407 }
1408 else {
1409 (*(g1++)) *= x * (*(g2++));
1410 if( ( k = *(++B2) ) == Inf< Index >() )
1411 break;
1412 if( ( h = *(++B1) ) == Inf< Index >() )
1413 break;
1414 }
1415 }
1416 } // end( VectScale( g1 , B1 , x , g2 , B2 )
1417
1418/*--------------------------------------------------------------------------*/
1419
1420template< class T1, class T2 >
1421inline void VectScaleB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1422{
1423 // g1{B} *= g2, all other entries of g1 unchanged
1424
1425 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1426 g1[ h ] *= *(g2++);
1427 }
1428
1429/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1430
1431template< class T1 , class T2 , class T3 >
1432inline void VectScaleB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1433{
1434 // g1{B} *= x * g2, all other entries of g1 unchanged
1435
1436 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1437 g1[ h ] *= x * (*(g2++));
1438 }
1439
1440/*--------------------------------------------------------------------------*/
1441
1442template< class T1, class T2 >
1443inline void VectScaleBB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1444{
1445 // g1{B} *= g2{B}, all other entries unchanged
1446
1447 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1448 g1[ h ] *= g2[ h ];
1449 }
1450
1451/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1452
1453template< class T1 , class T2 , class T3 >
1454inline void VectScaleBB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1455{
1456 // g1{B} *= x * g2{B}, all other entries unchanged
1457
1458 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1459 g1[ h ] *= x * g2[ h ];
1460 }
1461
1462/*--------------------------------------------------------------------------*/
1463/*--------------------------------------------------------------------------*/
1464
1465template< class T >
1466inline void VectIScale( T *g , const T x , Index n )
1467{
1468 // g := g / x, x a scalar; useful e.g. for *integer* types where
1469 // VectScale( g , 1 / x , n ) would not work
1470
1471 for( ; n-- ; )
1472 *(g++) /= x;
1473 }
1474
1475/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1476
1477template< class T >
1478inline void VectIScale( T *g , const T x , cIndex_Set B )
1479{
1480 // g{B} := g{B} / x, x a scalar
1481
1482 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1483 g[ h ] /= x;
1484 }
1485
1486/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1487
1488template< class T1, class T2 >
1489inline void VectIScale( T1 *g1 , const T2 *g2 , Index n )
1490{
1491 // g1 /= g2
1492
1493 for( ; n-- ; )
1494 *(g1++) /= *(g2++);
1495 }
1496
1497/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1498
1499template< class T1, class T2 >
1500inline void VectIScale( T1 *g1 , const T2 *g2 , cIndex_Set B )
1501{
1502 // g1 /= g2{B}
1503
1504 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1505 *(g1++) /= g2[ h ];
1506 }
1507
1508/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1509
1510template< class T1, class T2 >
1511inline void VectIScale( T1 *g1 , cIndex_Set B1 , const T2 *g2 , cIndex_Set B2 )
1512{
1513 // g1{B} /= g2{B}, B = intersection of B1 and B2, all other unchanged
1514
1515 Index h = *B1;
1516 Index k = *B2;
1517 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
1518 for(;;)
1519 if( k < h ) {
1520 if( ( k = *(++B2) ) == Inf< Index >() )
1521 break;
1522 g2++;
1523 }
1524 else
1525 if( h < k ) {
1526 if( ( h = *(++B1) ) == Inf< Index >() )
1527 break;
1528 g1++;
1529 }
1530 else {
1531 (*(g1++)) /= (*(g2++));
1532 if( ( k = *(++B2) ) == Inf< Index >() )
1533 break;
1534 if( ( h = *(++B1) ) == Inf< Index >() )
1535 break;
1536 }
1537 }
1538 } // end( VectIScale( g1 , B1 , g2 , B2 )
1539
1540/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1541
1542template< class T1 , class T2 , class T3 >
1543inline void VectIScale( T1 *g1 , const T2 *g2 , const T3 x , Index n )
1544{
1545 // g1 /= x * g2
1546
1547 for( ; n-- ; )
1548 *(g1++) /= ( x * (*(g2++)) );
1549 }
1550
1551/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1552
1553template< class T1 , class T2 , class T3 >
1554inline void VectIScale( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1555{
1556 // g1 /= x * g2{B}
1557
1558 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1559 *(g1++) /= ( x * g2[ h ] );
1560 }
1561
1562/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1563
1564template< class T1 , class T2 , class T3 >
1565inline void VectIScale( T1 *g1 , cIndex_Set B1 , const T3 x , const T2 *g2 ,
1566 cIndex_Set B2 )
1567{
1568 // g1{B} /= x * g2{B}, B = intersection of B1 and B2, all other unchanged
1569
1570 Index h = *B1;
1571 Index k = *B2;
1572 if( ( h < Inf< Index >() ) && ( k < Inf< Index >() ) ) {
1573 for(;;)
1574 if( k < h ) {
1575 if( ( k = *(++B2) ) == Inf< Index >() )
1576 break;
1577 g2++;
1578 }
1579 else
1580 if( h < k ) {
1581 if( ( h = *(++B1) ) == Inf< Index >() )
1582 break;
1583 g1++;
1584 }
1585 else {
1586 (*(g1++)) /= ( x * (*(g2++)) );
1587 if( ( k = *(++B2) ) == Inf< Index >() )
1588 break;
1589 if( ( h = *(++B1) ) == Inf< Index >() )
1590 break;
1591 }
1592 }
1593 } // end( VectIScale( g1 , B1 , x , g2 , B2 )
1594
1595/*--------------------------------------------------------------------------*/
1596
1597template< class T1, class T2 >
1598inline void VectIScaleB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1599{
1600 // g1{B} /= g2, all other entries of g1 unchanged
1601
1602 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1603 g1[ h ] /= *(g2++);
1604 }
1605
1606/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1607
1608template< class T1 , class T2 , class T3 >
1609inline void VectIScaleB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1610{
1611 // g1{B} /= x * g2, all other entries of g1 unchanged
1612
1613 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1614 g1[ h ] /= x * (*(g2++));
1615 }
1616
1617/*--------------------------------------------------------------------------*/
1618
1619template< class T1, class T2 >
1620inline void VectIScaleBB( T1 *g1 , const T2 *g2 , cIndex_Set B )
1621{
1622 // g1{B} /= g2{B}, all other entries unchanged
1623
1624 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1625 g1[ h ] *= g2[ h ];
1626 }
1627
1628/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1629
1630template< class T1 , class T2 , class T3 >
1631inline void VectIScaleBB( T1 *g1 , const T2 *g2 , const T3 x , cIndex_Set B )
1632{
1633 // g1{B} /= x * g2{B}, all other entries unchanged
1634
1635 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1636 g1[ h ] *= x * g2[ h ];
1637 }
1638
1639/*--------------------------------------------------------------------------*/
1640/*--------------------------------------------------------------------------*/
1641
1642template< class T1, class T2 >
1643inline void VectAdd( T1 *g , const T2 *g1 , const T1 x , Index n )
1644{
1645 // g := g1 + x
1646
1647 for( ; n-- ; )
1648 *(g++) = *(g1++) + x;
1649 }
1650
1651/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1652
1653template< class T1 , class T2 , class T3 >
1654inline void VectAdd( T1 *g , const T2 *g1 , const T3 *g2 , Index n )
1655{
1656 // g := g1 + g2
1657
1658 for( ; n-- ; )
1659 *(g++) = *(g1++) + *(g2++);
1660 }
1661
1662/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1663
1664template< class T1 , class T2 , class T3 >
1665inline void VectAdd( T1 *g , const T2 *g1 , const T3 *g2 , cIndex_Set B )
1666{
1667 // g := g1 + g2{B}
1668
1669 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1670 *(g++) = *(g1++) + g2[ h ];
1671 }
1672
1673/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1674
1675template< class T1 , class T2 , class T3 >
1676inline void VectAdd( T1 *g , const T2 *g1 , const T3 *g2 , const T3 x ,
1677 Index n )
1678{
1679 // g := g1 + x * g2
1680
1681 for( ; n-- ; )
1682 *(g++) = *(g1++) + x * (*(g2++));
1683 }
1684
1685/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1686
1687template< class T1, class T2 >
1688inline void VectAdd( T1 *g , const T2 *g1 , const T2 x1 , const T1 x , Index n )
1689{
1690 // g := x1 * g1 + x
1691
1692 for( ; n-- ; )
1693 *(g++) = x1 * (*(g1++)) + x;
1694 }
1695
1696/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1697
1698template< class T1 , class T2 , class T3 >
1699inline void VectAdd( T1 *g , const T2 *g1 , const T2 x1 , const T3 *g2 ,
1700 const T3 x2 , Index n )
1701{
1702 // g := x1 * g1 + x2 * g2
1703
1704 for( ; n-- ; )
1705 *(g++) = x1 * (*(g1++)) + x2 * (*(g2++));
1706 }
1707
1708/*--------------------------------------------------------------------------*/
1709/*--------------------------------------------------------------------------*/
1710
1711template< class T1 , class T2 , class T3 >
1712inline void VectDiff( T1 *g , const T2 *g1 , const T3 *g2 , Index n )
1713{
1714 // g := g1 - g2
1715
1716 for( ; n-- ; )
1717 *(g++) = *(g1++) - *(g2++);
1718 }
1719
1720/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1721
1722template< class T1 , class T2 , class T3 >
1723inline void VectDiff( T1 *g , const T2 *g1 , const T3 *g2 , cIndex_Set B )
1724{
1725 // g := g1 - g2{B}
1726
1727 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1728 *(g++) = *(g1++) - g2[ h ];
1729 }
1730
1731/*--------------------------------------------------------------------------*/
1732/*--------------------------------------------------------------------------*/
1733
1734template< class T1, class T2 >
1735inline void VectMult( T1 *g , const T2 *g1 , const T1 x , Index n )
1736{
1737 // g := g1 * x
1738
1739 for( ; n-- ; )
1740 *(g++) = *(g1++) * x;
1741 }
1742
1743/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1744
1745template< class T1 , class T2 , class T3 >
1746inline void VectMult( T1 *g , const T2 *g1 , const T3 *g2 , Index n )
1747{
1748 // g := g1 * g2
1749
1750 for( ; n-- ; )
1751 *(g++) = *(g1++) * (*(g2++));
1752 }
1753
1754/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1755
1756template< class T1 , class T2 , class T3 >
1757inline void VectMult( T1 *g , const T2 *g1 , const T3 *g2 , cIndex_Set B )
1758{
1759 // g := g1 * g2{B}
1760
1761 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1762 *(g++) = *(g1++) * g2[ h ];
1763 }
1764
1765/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1766
1767template< class T1 , class T2 , class T3 , class T4 >
1768inline void VectMult( T1 *g , const T2 x , const T3 *g1 , const T4 *g2 ,
1769 Index n )
1770{
1771 // g := x * g1 * g2
1772
1773 for( ; n-- ; )
1774 *(g++) = x * (*(g1++)) * (*(g2++));
1775 }
1776
1777/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1778
1779template< class T1 , class T2 , class T3 , class T4 >
1780inline void VectMult( T1 *g , const T2 x , const T3 *g1 , const T4 *g2 ,
1781 cIndex_Set B )
1782{
1783 // g := x * g1 * g2{B}
1784
1785 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1786 *(g++) = x * (*(g1++)) * g2[ h ];
1787 }
1788
1789/*--------------------------------------------------------------------------*/
1790/*--------------------------------------------------------------------------*/
1791
1792template< class T1 , class T2 , class T3 >
1793inline void VectDivide( T1 *g , const T2 *g1 , const T3 *g2 , Index n )
1794{
1795 // g := g1 / g2
1796
1797 for( ; n-- ; )
1798 *(g++) = *(g1++) / (*(g2++));
1799 }
1800
1801/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1802
1803template< class T1 , class T2 , class T3 >
1804inline void VectDivide( T1 *g , const T2 *g1 , const T3 *g2 , cIndex_Set B )
1805{
1806 // g := g1 / g2{B}
1807
1808 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1809 *(g++) = *(g1++) / g2[ h ];
1810 }
1811
1812/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1813
1814template< class T1 , class T2 , class T3 , class T4 >
1815inline void VectDivide( T1 *g , const T2 x , const T3 *g1 , const T4 *g2 ,
1816 Index n )
1817{
1818 // g := x * ( g1 / g2 )
1819
1820 for( ; n-- ; )
1821 *(g++) = x * ( (*(g1++)) / (*(g2++)) );
1822 }
1823
1824/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1825
1826template< class T1 , class T2 , class T3 , class T4 >
1827inline void VectDivide( T1 *g , const T2 x , const T3 *g1 , const T4 *g2 ,
1828 cIndex_Set B )
1829{
1830 // g := x * ( g1 / g2{B} )
1831
1832 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1833 *(g++) = x * ( (*(g1++)) / g2[ h ] );
1834 }
1835
1836/*--------------------------------------------------------------------------*/
1837/*--------------------------------------------------------------------------*/
1838
1839template< class T >
1840inline void VectXcg( T *g1 , T *g2 , Index n )
1841{
1842 // swap of g1 and g2 (element-wise)
1843
1844 for( ; n-- ; g1++ , g2++ )
1845 Swap( *g1 , *g2 );
1846 }
1847
1848/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1849
1850template< class T >
1851inline void VectXcg( T *g1 , T *g2 , cIndex_Set B )
1852{
1853 // swap of g1 and g2{B}
1854
1855 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; g1++ )
1856 Swap( *g1 , g2[ h ] );
1857 }
1858
1859/*--------------------------------------------------------------------------*/
1860
1861template< class T >
1862inline void VectXcgB( T *g1 , T *g2 , cIndex_Set B )
1863{
1864 // swap of g1{B} and g2
1865
1866 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; g2++ )
1867 Swap( g1[ h ] , *g2 );
1868 }
1869
1870/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
1871
1872template< class T >
1873inline void VectXcgBB( T *g1 , T *g2 , cIndex_Set B )
1874{
1875 // swap of g1{B} and g2{B}, all other entries are unchanged
1876
1877 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; )
1878 Swap( g1[ h ] , g2[ h ] );
1879 }
1880
1881/*--------------------------------------------------------------------------*/
1882/*--------------------------------------------------------------------------*/
1883/*-- Sparse/dense vector transformations: turn "sparse" vectors into --*/
1884/*-- "dense" ones and vice-versa. --*/
1885/*--------------------------------------------------------------------------*/
1886/*--------------------------------------------------------------------------*/
1887
1888template< class T >
1889inline Index_Set Sparsify( T* g , Index_Set B , Index n , Index Bs = 0 )
1890{
1891 // turns g from a "dense" n-vector to a "sparse" one, eliminating all items
1892 // that are exactly == 0; writes the set of nonzero items in B, with names
1893 // from Bs onwards, ordered in increasing sense; returns a pointer to the
1894 // first element in B after the last index written: this can be used for
1895 // computing the number of nonzeroes in the "sparsified" vector and/or for
1896 // Inf< Index >()-terminating the set
1897
1898 for( ; n ; n-- , g++ )
1899 if( *g )
1900 *(B++) = Bs++;
1901 else
1902 break;
1903
1904 if( n ) {
1905 T* tg = g++;
1906 for( Bs++ ; --n ; g++ , Bs++ )
1907 if( *g ) {
1908 *(tg++) = *g;
1909 *(B++) = Bs;
1910 }
1911 }
1912
1913 return( B );
1914
1915 } // end( Sparsify )
1916
1917/*--------------------------------------------------------------------------*/
1918
1919template< class T >
1920inline Index_Set SparsifyT( T* g , Index_Set B , Index n , const T eps ,
1921 Index Bs = 0 )
1922{
1923 // as Sparsify(), but elements are considered nonzero only if they are >=
1924 // eps (the idea is that all elements are >= 0)
1925
1926 for( ; n ; n-- , g++ )
1927 if( *g >= eps )
1928 *(B++) = Bs++;
1929 else
1930 break;
1931
1932 if( n ) {
1933 T* tg = g++;
1934 for( Bs++ ; --n ; g++ , Bs++ )
1935 if( *g >= eps ) {
1936 *(tg++) = *g;
1937 *(B++) = Bs;
1938 }
1939 }
1940
1941 return( B );
1942
1943 } // end( SparsifyT )
1944
1945/*--------------------------------------------------------------------------*/
1946
1947template< class T >
1948inline Index_Set SparsifyAT( T* g , Index_Set B , Index n , const T eps ,
1949 Index Bs = 0 )
1950{
1951 // as SparsifyT(), but elements are considered nonzero only if their ABS()
1952 // is >= eps
1953
1954 for( ; n ; n-- , g++ )
1955 if( ABS( *g ) >= eps )
1956 *(B++) = Bs++;
1957 else
1958 break;
1959
1960 if( n ) {
1961 T* tg = g++;
1962 for( Bs++ ; --n ; g++ , Bs++ )
1963 if( ABS( *g ) >= eps ) {
1964 *(tg++) = *g;
1965 *(B++) = Bs;
1966 }
1967 }
1968
1969 return( B );
1970
1971 } // end( SparsifyAT )
1972
1973/*--------------------------------------------------------------------------*/
1974
1975template< class T >
1976inline void Densify( T* g , cIndex_Set B , Index m , Index n , cIndex k = 0 )
1977{
1978 // turns g from a "sparse" m-vector, whose set of nonzero elements is B, to
1979 // a "dense" n-vector padded with zeroes where necessary; B has to be ordered
1980 // in increasing sense, but does not need to be Inf< Index >()-terminated (m
1981 // gives the same information); note that the function will write in
1982 // g[ n - 1 ], hence the vector has to have been properly allocated
1983 // if k > 0, the function only works for the subvector of g between k and
1984 // n - 1, that is, g[ 0 ] .. g[ k - 1 ] are left intact, while the
1985 // subvector is densified; it is required that B[] only contains indices
1986 // >= k (and < n, of course)
1987
1988 // since the vector is densified in place, processing has to be done from
1989 // the end backwards to avoid overriding the data
1990 // note that there are two different cases:
1991 // - if B[ 0 ] > k, all the elements g[ k ] .. g[ B[ 0 ] - 1 ] are zero, and
1992 // the beginning of B[] will be hit before hitting the g + k
1993 // - if instead B[ 0 ] == k, then for some h >= 0 we have that
1994 // g[ k ] .. g[ k + h ] are "already in the right place" (the beginning
1995 // of the "sparse" g[] is in fact dense already), and therefore one can
1996 // stop copying at g + h, before hitting the beginning of B[]
1997 // we distinguish between the two cases in order to optimize the code
1998
1999 if( m ) { // there is at least a nonzero element
2000 T *g2 = g + n;
2001 T *g1 = g + k + m - 1;
2002
2003 if( *B > k ) { // first case - - - - - - - - - - - - - - - - - - - - - - -
2004
2005 cIndex_Set B1 = B + m;
2006 for( Index h = *(--B1) ;; )
2007 if( h == --g2 - g ) {
2008 *g2 = *(g1--);
2009 if( B1 > B )
2010 h = *(--B1);
2011 else
2012 break;
2013 }
2014 else
2015 *g2 = 0;
2016
2017 for( ; --g2 > g1 ; )
2018 *g2 = 0;
2019 }
2020 else { // second case- - - - - - - - - - - - - - - - - - - - - - -
2021
2022 B += m;
2023 for( Index h = *(--B) ; --g2 > g1 ; )
2024 if( h == g2 - g ) {
2025 *g2 = *(g1--);
2026 h = *(--B);
2027 }
2028 else
2029 *g2 = 0;
2030
2031 } // end second case- - - - - - - - - - - - - - - - - - - - -
2032 }
2033 else // m == 0: just zero-out the vector
2034 VectAssign( g + k , T( 0 ) , n - k );
2035
2036 } // end( Densify )
2037
2038/*--------------------------------------------------------------------------*/
2039
2040template< class T >
2041inline void Compact( T* g , cIndex_Set B , Index n )
2042{
2043 // takes a "dense" n-vector g and "compacts" it deleting the elements whose
2044 // indices are in B (all elements of B[] must be in the range 0 .. n, B[]
2045 // must be ordered in increasing sense and Inf< Index >()-terminated)
2046 // the remaining entries in g[] are shifted left of the minimum possible
2047 // amount in order to fill the holes left by the deleted ones
2048
2049 Index i = *(B++); // current position where to write
2050 Index j = i + 1; // element to copy
2051
2052 for( Index h = *(B++) ; h < Inf< Index >() ; j++ ) {
2053 while( j < h )
2054 g[ i++ ] = g[ j++ ];
2055
2056 h = *(B++);
2057 }
2058
2059 VectAssign( g + i , g + j , n - j );
2060
2061 } // end( Compact )
2062
2063/*--------------------------------------------------------------------------*/
2064/*--------------------------------------------------------------------------*/
2065/*-- --*/
2066/*-- Functions for array manipulation (merging, shifting ...). The types --*/
2067/*-- involved in these operations must support ordering relations '==' --*/
2068/*-- and '>', as well as the assigment '='. --*/
2069/*-- --*/
2070/*--------------------------------------------------------------------------*/
2071/*--------------------------------------------------------------------------*/
2072
2073template< class T >
2074inline void Merge( T *g , const T *g2 , const T *g3 , const T Stp )
2075{
2076 // merges the two ordered and "Stp-terminated" vectors g2 and g3 (that is,
2077 // there must be an element in both g2 and g3 that is >= Stp, and that is
2078 // is considered to be the terminator) into g1, that will therefore be
2079 // ordered and Stp-terminated
2080
2081 T i = *g2;
2082 T j = *g3;
2083
2084 for(;;)
2085 if( i < j ) {
2086 if( ( *(g++) = i ) >= Stp )
2087 break;
2088
2089 i = *(++g2);
2090 }
2091 else {
2092 if( ( *(g++) = j ) >= Stp )
2093 break;
2094
2095 if( i == j )
2096 i = *(++g2);
2097
2098 j = *(++g3);
2099 }
2100 } // end( Merge )
2101
2102/*--------------------------------------------------------------------------*/
2103/*--------------------------------------------------------------------------*/
2104
2105template< class T >
2106inline void ShiftVect( T *g , Index n )
2107{
2108 // g[ 0 ] = g[ 1 ], g[ 1 ] = g[ 2 ] ... g[ n - 1 ] = g[ n ]
2109
2110 for( T *t = g ; n-- ; t = g )
2111 *t = *(++g);
2112 }
2113
2114/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
2115
2116template< class T >
2117inline void ShiftVect( T *g , Index n , cIndex k )
2118{
2119 // g[ i ] = g[ i + k ], i = 0 ... n - 1
2120
2121 for( T *t = g + k ; n-- ; )
2122 *(g++) = *(t++);
2123 }
2124
2125/*--------------------------------------------------------------------------*/
2126
2127template< class T >
2128inline void RotateVect( T *g , Index n )
2129{
2130 // as ShiftVect(), but g[ 0 ] is moved to g[ n ]
2131
2132 const T tmp = *g;
2133
2134 for( T *t = g ; n-- ; t = g )
2135 *t = *(++g);
2136
2137 *g = tmp;
2138 }
2139
2140/*--------------------------------------------------------------------------*/
2141
2142template< class T >
2143inline void ShiftRVect( T *g , Index n )
2144{
2145 // g[ n ] = g[ n - 1 ], g[ n - 1 ] = g[ n - 2 ], ..., g[ 1 ] = g[ 0 ]
2146
2147 for( T *t = ( g += n ) ; n-- ; t = g )
2148 *t = *(--g);
2149 }
2150
2151/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
2152
2153template< class T >
2154inline void ShiftRVect( T *g , Index n , cIndex k )
2155{
2156 // g[ i ] = g[ i - k ], i = n - 1 ... 0
2157
2158 g += n;
2159 for( T *t = g + k ; n-- ; )
2160 *(--t) = *(--g);
2161 }
2162
2163/*--------------------------------------------------------------------------*/
2164
2165template< class T >
2166inline void RotateRVect( T *g , Index n )
2167{
2168 // as ShiftRVect(), but g[ n ] is moved to g[ 0 ]
2169
2170 T *t = ( g += n );
2171 const T tmp = *t;
2172
2173 for( ; n-- ; t = g )
2174 *t = *(--g);
2175
2176 *g = tmp;
2177 }
2178
2179/*--------------------------------------------------------------------------*/
2180/*--------------------------------------------------------------------------*/
2181/*-- --*/
2182/*-- Searching/ordering functions. The types involved in these operations --*/
2183/*-- must support ordering relations '==' and '>', as well as the --*/
2184/*-- assigment '='. --*/
2185/*-- --*/
2186/*--------------------------------------------------------------------------*/
2187/*--------------------------------------------------------------------------*/
2188
2189template< class T >
2190inline Index Match( const T *g , const T x , cIndex n )
2191{
2192 // if v contains elements identical to x, then the smallest index among all
2193 // such elements is reported (0 .. n - 1), else n is reported
2194
2195 Index i = 0;
2196 for( ; ( i < n ) && ( *(g++) != x ) ; )
2197 i++;
2198
2199 return( i );
2200 }
2201
2202/*--------------------------------------------------------------------------*/
2203
2204template< class T >
2205inline bool EqualVect( const T *g1 , const T *g2 , Index n )
2206{
2207 // returns true <=> the two n-vectors g1 and g2 are element-wise identical
2208
2209 for( ; n-- ; )
2210 if( *(g1++) != *(g2++) )
2211 return( false );
2212
2213 return( true );
2214 }
2215
2216/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
2217
2218template< class T >
2219inline bool EqualVect( const T *g1 , const T *g2 , cIndex n , cIndex_Set B )
2220{
2221 // returns true <=> the two n-vectors g1 and g2 are element-wise identical,
2222 // where g2 is given in sparse form, i.e., g2[ i ] is the B[ i ]-th element
2223
2224 Index i = 0;
2225 for( Index h ; ( h = *(B++) ) < Inf< Index >() ; i++ ) {
2226 for( ; i < h ; i++ )
2227 if( *(g1++) )
2228 return( false );
2229
2230 if( *(g1++) != *(g2++) )
2231 return( false );
2232 }
2233
2234 for( ; i < n ; i++ )
2235 if( *(g1++) )
2236 return( false );
2237
2238 return( true );
2239 }
2240
2241/*--------------------------------------------------------------------------*/
2242/*--------------------------------------------------------------------------*/
2243
2244template< class T >
2245inline Index BinSearch( const T *Set , Index Stop , const T What )
2246{
2247 // perform a binary search on the ordered set of Stop elements (of type T and
2248 // without replications) contained in the array Set: Set must be "infinity-
2249 // terminated", i.e. Set[ Stop ] > What. Searches for the element What, that
2250 // may or may not be in Set: if What is found, its position is reported,
2251 // otherwise the position of the smallest element > What in Set is reported.
2252 // If there are no elements > What in Set, then Stop is reported
2253
2254 Index i = Stop / 2;
2255 for( Index Strt = 0 ; Strt != Stop ; ) {
2256 if( Set[ i ] == What )
2257 break;
2258 else
2259 if( Set[ i ] > What )
2260 Stop = i;
2261 else
2262 Strt = i + 1;
2263
2264 i = ( Strt + Stop ) / 2;
2265 }
2266
2267 return( i );
2268
2269 } // end( BinSearch )
2270
2271/*--------------------------------------------------------------------------*/
2272
2273template< class T >
2274inline Index BinSearch1( const T *Set , Index Stop , const T What )
2275{
2276 // like BinSearch() above, but What *must be* in Set (=> Set is nonempty)
2277
2278 Index Strt = 0;
2279 Index i = Stop / 2;
2280 for( Index h ; ( h = Set[ i ] ) != What ; ) {
2281 if( h > What )
2282 Stop = i - 1;
2283 else
2284 Strt = i + 1;
2285
2286 i = ( Strt + Stop ) / 2;
2287 }
2288
2289 return( i );
2290
2291 } // end( BinSearch1 )
2292
2293/*--------------------------------------------------------------------------*/
2294
2295template< class T >
2296inline Index BinSearch2( const T *Set , Index Stop , const T What )
2297{
2298 // like BinSearch() above, but What must *not* be in Set (=> the position
2299 // of the smallest element > What in Set will be returned)
2300
2301 Index i = Stop / 2;
2302 for( Index Strt = 0 ; Strt != Stop ; ) {
2303 if( Set[ i ] > What )
2304 Stop = i;
2305 else
2306 Strt = i + 1;
2307
2308 i = ( Strt + Stop ) / 2;
2309 }
2310
2311 return( i );
2312
2313 } // end( BinSearch2 )
2314
2315/*--------------------------------------------------------------------------*/
2316/*--------------------------------------------------------------------------*/
2317
2318template< class T >
2319inline void HeapIns( T *H , const T x , Index n )
2320{
2321 // H is a binary heap of elements of type T, ordered in increasing sense
2322 // (i.e. the root of the heap is the smallest element) and containing n
2323 // elements: the new element x is inserted in H
2324
2325 for( H-- , n++ ;;) {
2326 Index p = n / 2;
2327 if( ( ! p ) || ( H[ p ] <= x ) )
2328 break;
2329
2330 H[ n ] = H[ p ];
2331 n = p;
2332 }
2333
2334 H[ n ] = x;
2335
2336 } // end( HeapIns )
2337
2338/*--------------------------------------------------------------------------*/
2339
2340template< class T >
2341inline Index HeapDel( T *H , cIndex n )
2342{
2343 // H is a binary heap of elements of type T, as above: returns the smallest
2344 // element (the root) deleting it from H. n is now intended to be the
2345 // position of the last element in H, i.e. | H | - 1, i.e. | H | *after*
2346 // the deletion
2347
2348 const T h = *H;
2349 *H = H[ n ];
2350
2351 for( Index i = 0 ; i < n ; ) {
2352 Index j = 2 * i;
2353
2354 if( ++j >= n )
2355 break;
2356
2357 Index k = j;
2358
2359 if( ++k < n )
2360 if( H[ k ] < H[ j ] )
2361 k = j;
2362
2363 if( H[ i ] <= H[ j ] )
2364 break;
2365
2366 Swap( H[ i ] , H[ j ] );
2367 i = j;
2368 }
2369
2370 return( h );
2371 }
2372
2373/*--------------------------------------------------------------------------*/
2374
2375}; // end( namespace OPTtypes_di_unipi_it )
2376
2377/*--------------------------------------------------------------------------*/
2378/*--------------------------------------------------------------------------*/
2379
2380#endif /* OPTvect.h included */
2381
2382/*--------------------------------------------------------------------------*/
2383/*---------------------- End File OPTvect.h --------------------------------*/
2384/*--------------------------------------------------------------------------*/
OPTtypes_di_unipi_it::cIndex
const Index cIndex
a read-only Index
Definition OPTtypes.h:64
OPTtypes_di_unipi_it::cIndex_Set
cIndex * cIndex_Set
read-only array
Definition OPTtypes.h:65
OPTtypes_di_unipi_it::Index_Set
Index * Index_Set
set (array) of indices
Definition OPTtypes.h:61
OPTtypes_di_unipi_it::Index
unsigned int Index
Index in a vector ( >= 0 )
Definition OPTtypes.h:60
OPTtypes_di_unipi_it::Inf
static constexpr T Inf(void) noexcept
Inf< T >() = infinity value for T.
Definition OPTUtils.h:357
OPTtypes_di_unipi_it
Definition OPTtypes.h:41
Generated by doxygen 1.12.0