c890253ec8fb5461ce51ba93e29f2e4a384d03e8
[scilab.git] / scilab / modules / scicos / src / scicos_sundials / src / nvec_ser / nvector_serial.c
1 /*
2  * -----------------------------------------------------------------
3  * $Revision: 1.1 $
4  * $Date: 2006/07/05 15:32:37 $
5  * ----------------------------------------------------------------- 
6  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, Radu Serban,
7  *                and Aaron Collier @ LLNL
8  * -----------------------------------------------------------------
9  * Copyright (c) 2002, The Regents of the University of California.
10  * Produced at the Lawrence Livermore National Laboratory.
11  * All rights reserved.
12  * For details, see the LICENSE file.
13  * -----------------------------------------------------------------
14  * This is the implementation file for a serial implementation
15  * of the NVECTOR package.
16  * -----------------------------------------------------------------
17  */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21
22 #include <nvector/nvector_serial.h>
23 #include <sundials/sundials_math.h>
24
25 #define ZERO   RCONST(0.0)
26 #define HALF   RCONST(0.5)
27 #define ONE    RCONST(1.0)
28 #define ONEPT5 RCONST(1.5)
29
30 /* Private function prototypes */
31 /* z=x */
32 static void VCopy_Serial(N_Vector x, N_Vector z);
33 /* z=x+y */
34 static void VSum_Serial(N_Vector x, N_Vector y, N_Vector z);
35 /* z=x-y */
36 static void VDiff_Serial(N_Vector x, N_Vector y, N_Vector z);
37 /* z=-x */
38 static void VNeg_Serial(N_Vector x, N_Vector z);
39 /* z=c(x+y) */
40 static void VScaleSum_Serial(realtype c, N_Vector x, N_Vector y, N_Vector z);
41 /* z=c(x-y) */
42 static void VScaleDiff_Serial(realtype c, N_Vector x, N_Vector y, N_Vector z); 
43 /* z=ax+y */
44 static void VLin1_Serial(realtype a, N_Vector x, N_Vector y, N_Vector z);
45 /* z=ax-y */
46 static void VLin2_Serial(realtype a, N_Vector x, N_Vector y, N_Vector z);
47 /* y <- ax+y */
48 static void Vaxpy_Serial(realtype a, N_Vector x, N_Vector y);
49 /* x <- ax */
50 static void VScaleBy_Serial(realtype a, N_Vector x);
51
52 /*
53  * -----------------------------------------------------------------
54  * exported functions
55  * -----------------------------------------------------------------
56  */
57
58 /* ----------------------------------------------------------------------------
59  * Function to create a new empty serial vector 
60  */
61
62 N_Vector N_VNewEmpty_Serial(long int length)
63 {
64   N_Vector v;
65   N_Vector_Ops ops;
66   N_VectorContent_Serial content;
67
68   /* Create vector */
69   v = NULL;
70   v = (N_Vector) malloc(sizeof *v);
71   if (v == NULL) return(NULL);
72   
73   /* Create vector operation structure */
74   ops = NULL;
75   ops = (N_Vector_Ops) malloc(sizeof(struct _generic_N_Vector_Ops));
76   if (ops == NULL) { free(v); return(NULL); }
77
78   ops->nvclone           = N_VClone_Serial;
79   ops->nvcloneempty      = N_VCloneEmpty_Serial;
80   ops->nvdestroy         = N_VDestroy_Serial;
81   ops->nvspace           = N_VSpace_Serial;
82   ops->nvgetarraypointer = N_VGetArrayPointer_Serial;
83   ops->nvsetarraypointer = N_VSetArrayPointer_Serial;
84   ops->nvlinearsum       = N_VLinearSum_Serial;
85   ops->nvconst           = N_VConst_Serial;
86   ops->nvprod            = N_VProd_Serial;
87   ops->nvdiv             = N_VDiv_Serial;
88   ops->nvscale           = N_VScale_Serial;
89   ops->nvabs             = N_VAbs_Serial;
90   ops->nvinv             = N_VInv_Serial;
91   ops->nvaddconst        = N_VAddConst_Serial;
92   ops->nvdotprod         = N_VDotProd_Serial;
93   ops->nvmaxnorm         = N_VMaxNorm_Serial;
94   ops->nvwrmsnormmask    = N_VWrmsNormMask_Serial;
95   ops->nvwrmsnorm        = N_VWrmsNorm_Serial;
96   ops->nvmin             = N_VMin_Serial;
97   ops->nvwl2norm         = N_VWL2Norm_Serial;
98   ops->nvl1norm          = N_VL1Norm_Serial;
99   ops->nvcompare         = N_VCompare_Serial;
100   ops->nvinvtest         = N_VInvTest_Serial;
101   ops->nvconstrmask      = N_VConstrMask_Serial;
102   ops->nvminquotient     = N_VMinQuotient_Serial;
103
104   /* Create content */
105   content = NULL;
106   content = (N_VectorContent_Serial) malloc(sizeof(struct _N_VectorContent_Serial));
107   if (content == NULL) { free(ops); free(v); return(NULL); }
108
109   content->length   = length;
110   content->own_data = FALSE;
111   content->data     = NULL;
112
113   /* Attach content and ops */
114   v->content = content;
115   v->ops     = ops;
116
117   return(v);
118 }
119
120 /* ----------------------------------------------------------------------------
121  * Function to create a new serial vector 
122  */
123
124 N_Vector N_VNew_Serial(long int length)
125 {
126   N_Vector v;
127   realtype *data;
128
129   v = NULL;
130   v = N_VNewEmpty_Serial(length);
131   if (v == NULL) return(NULL);
132
133   /* Create data */
134   if (length > 0) {
135
136     /* Allocate memory */
137     data = NULL;
138     data = (realtype *) malloc(length * sizeof(realtype));
139     if(data == NULL) { N_VDestroy_Serial(v); return(NULL); }
140
141     /* Attach data */
142     NV_OWN_DATA_S(v) = TRUE;
143     NV_DATA_S(v)     = data;
144
145   }
146
147   return(v);
148 }
149
150 /* ----------------------------------------------------------------------------
151  * Function to create a serial N_Vector with user data component 
152  */
153
154 N_Vector N_VMake_Serial(long int length, realtype *v_data)
155 {
156   N_Vector v;
157
158   v = NULL;
159   v = N_VNewEmpty_Serial(length);
160   if (v == NULL) return(NULL);
161
162   if (length > 0) {
163     /* Attach data */
164     NV_OWN_DATA_S(v) = FALSE;
165     NV_DATA_S(v)     = v_data;
166   }
167
168   return(v);
169 }
170
171 /* ----------------------------------------------------------------------------
172  * Function to create an array of new serial vectors. 
173  */
174
175 N_Vector *N_VCloneVectorArray_Serial(int count, N_Vector w)
176 {
177   N_Vector *vs;
178   int j;
179
180   if (count <= 0) return(NULL);
181
182   vs = NULL;
183   vs = (N_Vector *) malloc(count * sizeof(N_Vector));
184   if(vs == NULL) return(NULL);
185
186   for (j = 0; j < count; j++) {
187     vs[j] = NULL;
188     vs[j] = N_VClone_Serial(w);
189     if (vs[j] == NULL) {
190       N_VDestroyVectorArray_Serial(vs, j-1);
191       return(NULL);
192     }
193   }
194
195   return(vs);
196 }
197
198 /* ----------------------------------------------------------------------------
199  * Function to create an array of new serial vectors with NULL data array. 
200  */
201
202 N_Vector *N_VCloneVectorArrayEmpty_Serial(int count, N_Vector w)
203 {
204   N_Vector *vs;
205   int j;
206
207   if (count <= 0) return(NULL);
208
209   vs = NULL;
210   vs = (N_Vector *) malloc(count * sizeof(N_Vector));
211   if(vs == NULL) return(NULL);
212
213   for (j = 0; j < count; j++) {
214     vs[j] = NULL;
215     vs[j] = N_VCloneEmpty_Serial(w);
216     if (vs[j] == NULL) {
217       N_VDestroyVectorArray_Serial(vs, j-1);
218       return(NULL);
219     }
220   }
221
222   return(vs);
223 }
224
225 /* ----------------------------------------------------------------------------
226  * Function to free an array created with N_VCloneVectorArray_Serial
227  */
228
229 void N_VDestroyVectorArray_Serial(N_Vector *vs, int count)
230 {
231   int j;
232
233   for (j = 0; j < count; j++) N_VDestroy_Serial(vs[j]);
234
235   free(vs); vs = NULL;
236
237   return;
238 }
239
240 /* ----------------------------------------------------------------------------
241  * Function to print the a serial vector 
242  */
243  
244 void N_VPrint_Serial(N_Vector x)
245 {
246   long int i, N;
247   realtype *xd;
248
249   xd = NULL;
250
251   N  = NV_LENGTH_S(x);
252   xd = NV_DATA_S(x);
253
254   for (i = 0; i < N; i++) {
255 #if defined(SUNDIALS_EXTENDED_PRECISION)
256     printf("%11.8Lg\n", xd[i]);
257 #elif defined(SUNDIALS_DOUBLE_PRECISION)
258     printf("%11.8lg\n", xd[i]);
259 #else
260     printf("%11.8g\n", xd[i]);
261 #endif
262   }
263   printf("\n");
264
265   return;
266 }
267
268 /*
269  * -----------------------------------------------------------------
270  * implementation of vector operations
271  * -----------------------------------------------------------------
272  */
273
274 N_Vector N_VCloneEmpty_Serial(N_Vector w)
275 {
276   N_Vector v;
277   N_Vector_Ops ops;
278   N_VectorContent_Serial content;
279
280   if (w == NULL) return(NULL);
281
282   /* Create vector */
283   v = NULL;
284   v = (N_Vector) malloc(sizeof *v);
285   if (v == NULL) return(NULL);
286
287   /* Create vector operation structure */
288   ops = NULL;
289   ops = (N_Vector_Ops) malloc(sizeof(struct _generic_N_Vector_Ops));
290   if (ops == NULL) { free(v); return(NULL); }
291   
292   ops->nvclone           = w->ops->nvclone;
293   ops->nvcloneempty      = w->ops->nvcloneempty;
294   ops->nvdestroy         = w->ops->nvdestroy;
295   ops->nvspace           = w->ops->nvspace;
296   ops->nvgetarraypointer = w->ops->nvgetarraypointer;
297   ops->nvsetarraypointer = w->ops->nvsetarraypointer;
298   ops->nvlinearsum       = w->ops->nvlinearsum;
299   ops->nvconst           = w->ops->nvconst;  
300   ops->nvprod            = w->ops->nvprod;   
301   ops->nvdiv             = w->ops->nvdiv;
302   ops->nvscale           = w->ops->nvscale; 
303   ops->nvabs             = w->ops->nvabs;
304   ops->nvinv             = w->ops->nvinv;
305   ops->nvaddconst        = w->ops->nvaddconst;
306   ops->nvdotprod         = w->ops->nvdotprod;
307   ops->nvmaxnorm         = w->ops->nvmaxnorm;
308   ops->nvwrmsnormmask    = w->ops->nvwrmsnormmask;
309   ops->nvwrmsnorm        = w->ops->nvwrmsnorm;
310   ops->nvmin             = w->ops->nvmin;
311   ops->nvwl2norm         = w->ops->nvwl2norm;
312   ops->nvl1norm          = w->ops->nvl1norm;
313   ops->nvcompare         = w->ops->nvcompare;    
314   ops->nvinvtest         = w->ops->nvinvtest;
315   ops->nvconstrmask      = w->ops->nvconstrmask;
316   ops->nvminquotient     = w->ops->nvminquotient;
317
318   /* Create content */
319   content = NULL;
320   content = (N_VectorContent_Serial) malloc(sizeof(struct _N_VectorContent_Serial));
321   if (content == NULL) { free(ops); free(v); return(NULL); }
322
323   content->length   = NV_LENGTH_S(w);
324   content->own_data = FALSE;
325   content->data     = NULL;
326
327   /* Attach content and ops */
328   v->content = content;
329   v->ops     = ops;
330
331   return(v);
332 }
333
334 N_Vector N_VClone_Serial(N_Vector w)
335 {
336   N_Vector v;
337   realtype *data;
338   long int length;
339
340   v = NULL;
341   v = N_VCloneEmpty_Serial(w);
342   if (v == NULL) return(NULL);
343
344   length = NV_LENGTH_S(w);
345
346   /* Create data */
347   if (length > 0) {
348
349     /* Allocate memory */
350     data = NULL;
351     data = (realtype *) malloc(length * sizeof(realtype));
352     if(data == NULL) { N_VDestroy_Serial(v); return(NULL); }
353
354     /* Attach data */
355     NV_OWN_DATA_S(v) = TRUE;
356     NV_DATA_S(v)     = data;
357
358   }
359
360   return(v);
361 }
362
363 void N_VDestroy_Serial(N_Vector v)
364 {
365   if (NV_OWN_DATA_S(v) == TRUE) {
366     free(NV_DATA_S(v));
367     NV_DATA_S(v) = NULL;
368   }
369   free(v->content); v->content = NULL;
370   free(v->ops); v->ops = NULL;
371   free(v); v = NULL;
372
373   return;
374 }
375
376 void N_VSpace_Serial(N_Vector v, long int *lrw, long int *liw)
377 {
378   *lrw = NV_LENGTH_S(v);
379   *liw = 1;
380
381   return;
382 }
383
384 realtype *N_VGetArrayPointer_Serial(N_Vector v)
385 {
386   return((realtype *) NV_DATA_S(v));
387 }
388
389 void N_VSetArrayPointer_Serial(realtype *v_data, N_Vector v)
390 {
391   if (NV_LENGTH_S(v) > 0) NV_DATA_S(v) = v_data;
392
393   return;
394 }
395
396 void N_VLinearSum_Serial(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z)
397 {
398   long int i, N;
399   realtype c, *xd, *yd, *zd;
400   N_Vector v1, v2;
401   booleantype test;
402
403   xd = yd = zd = NULL;
404
405   if ((b == ONE) && (z == y)) {    /* BLAS usage: axpy y <- ax+y */
406     Vaxpy_Serial(a,x,y);
407     return;
408   }
409
410   if ((a == ONE) && (z == x)) {    /* BLAS usage: axpy x <- by+x */
411     Vaxpy_Serial(b,y,x);
412     return;
413   }
414
415   /* Case: a == b == 1.0 */
416
417   if ((a == ONE) && (b == ONE)) {
418     VSum_Serial(x, y, z);
419     return;
420   }
421
422   /* Cases: (1) a == 1.0, b = -1.0, (2) a == -1.0, b == 1.0 */
423
424   if ((test = ((a == ONE) && (b == -ONE))) || ((a == -ONE) && (b == ONE))) {
425     v1 = test ? y : x;
426     v2 = test ? x : y;
427     VDiff_Serial(v2, v1, z);
428     return;
429   }
430
431   /* Cases: (1) a == 1.0, b == other or 0.0, (2) a == other or 0.0, b == 1.0 */
432   /* if a or b is 0.0, then user should have called N_VScale */
433
434   if ((test = (a == ONE)) || (b == ONE)) {
435     c  = test ? b : a;
436     v1 = test ? y : x;
437     v2 = test ? x : y;
438     VLin1_Serial(c, v1, v2, z);
439     return;
440   }
441
442   /* Cases: (1) a == -1.0, b != 1.0, (2) a != 1.0, b == -1.0 */
443
444   if ((test = (a == -ONE)) || (b == -ONE)) {
445     c = test ? b : a;
446     v1 = test ? y : x;
447     v2 = test ? x : y;
448     VLin2_Serial(c, v1, v2, z);
449     return;
450   }
451
452   /* Case: a == b */
453   /* catches case both a and b are 0.0 - user should have called N_VConst */
454
455   if (a == b) {
456     VScaleSum_Serial(a, x, y, z);
457     return;
458   }
459
460   /* Case: a == -b */
461
462   if (a == -b) {
463     VScaleDiff_Serial(a, x, y, z);
464     return;
465   }
466
467   /* Do all cases not handled above:
468      (1) a == other, b == 0.0 - user should have called N_VScale
469      (2) a == 0.0, b == other - user should have called N_VScale
470      (3) a,b == other, a !=b, a != -b */
471   
472   N  = NV_LENGTH_S(x);
473   xd = NV_DATA_S(x);
474   yd = NV_DATA_S(y);
475   zd = NV_DATA_S(z);
476
477   for (i = 0; i < N; i++)
478     zd[i] = (a*xd[i])+(b*yd[i]);
479
480   return;
481 }
482
483 void N_VConst_Serial(realtype c, N_Vector z)
484 {
485   long int i, N;
486   realtype *zd;
487
488   zd = NULL;
489
490   N  = NV_LENGTH_S(z);
491   zd = NV_DATA_S(z);
492
493   for (i = 0; i < N; i++) zd[i] = c;
494
495   return;
496 }
497
498 void N_VProd_Serial(N_Vector x, N_Vector y, N_Vector z)
499 {
500   long int i, N;
501   realtype *xd, *yd, *zd;
502
503   xd = yd = zd = NULL;
504
505   N  = NV_LENGTH_S(x);
506   xd = NV_DATA_S(x);
507   yd = NV_DATA_S(y);
508   zd = NV_DATA_S(z);
509
510   for (i = 0; i < N; i++)
511     zd[i] = xd[i]*yd[i];
512
513   return;
514 }
515
516 void N_VDiv_Serial(N_Vector x, N_Vector y, N_Vector z)
517 {
518   long int i, N;
519   realtype *xd, *yd, *zd;
520
521   xd = yd = zd = NULL;
522
523   N  = NV_LENGTH_S(x);
524   xd = NV_DATA_S(x);
525   yd = NV_DATA_S(y);
526   zd = NV_DATA_S(z);
527
528   for (i = 0; i < N; i++)
529     zd[i] = xd[i]/yd[i];
530
531   return;
532 }
533
534 void N_VScale_Serial(realtype c, N_Vector x, N_Vector z)
535 {
536   long int i, N;
537   realtype *xd, *zd;
538
539   xd = zd = NULL;
540
541   if (z == x) {  /* BLAS usage: scale x <- cx */
542     VScaleBy_Serial(c, x);
543     return;
544   }
545
546   if (c == ONE) {
547     VCopy_Serial(x, z);
548   } else if (c == -ONE) {
549     VNeg_Serial(x, z);
550   } else {
551     N  = NV_LENGTH_S(x);
552     xd = NV_DATA_S(x);
553     zd = NV_DATA_S(z);
554     for (i = 0; i < N; i++) 
555       zd[i] = c*xd[i];
556   }
557
558   return;
559 }
560
561 void N_VAbs_Serial(N_Vector x, N_Vector z)
562 {
563   long int i, N;
564   realtype *xd, *zd;
565
566   xd = zd = NULL;
567
568   N  = NV_LENGTH_S(x);
569   xd = NV_DATA_S(x);
570   zd = NV_DATA_S(z);
571
572   for (i = 0; i < N; i++)
573     zd[i] = ABS(xd[i]);
574
575   return;
576 }
577
578 void N_VInv_Serial(N_Vector x, N_Vector z)
579 {
580   long int i, N;
581   realtype *xd, *zd;
582
583   xd = zd = NULL;
584
585   N  = NV_LENGTH_S(x);
586   xd = NV_DATA_S(x);
587   zd = NV_DATA_S(z);
588
589   for (i = 0; i < N; i++)
590     zd[i] = ONE/xd[i];
591
592   return;
593 }
594
595 void N_VAddConst_Serial(N_Vector x, realtype b, N_Vector z)
596 {
597   long int i, N;
598   realtype *xd, *zd;
599
600   xd = zd = NULL;
601
602   N  = NV_LENGTH_S(x);
603   xd = NV_DATA_S(x);
604   zd = NV_DATA_S(z);
605
606   for (i = 0; i < N; i++) 
607     zd[i] = xd[i]+b;
608
609   return;
610 }
611
612 realtype N_VDotProd_Serial(N_Vector x, N_Vector y)
613 {
614   long int i, N;
615   realtype sum, *xd, *yd;
616
617   sum = ZERO;
618   xd = yd = NULL;
619
620   N  = NV_LENGTH_S(x);
621   xd = NV_DATA_S(x);
622   yd = NV_DATA_S(y);
623
624   for (i = 0; i < N; i++)
625     sum += xd[i]*yd[i];
626   
627   return(sum);
628 }
629
630 realtype N_VMaxNorm_Serial(N_Vector x)
631 {
632   long int i, N;
633   realtype max, *xd;
634
635   max = ZERO;
636   xd = NULL;
637
638   N  = NV_LENGTH_S(x);
639   xd = NV_DATA_S(x);
640
641   for (i = 0; i < N; i++) {
642     if (ABS(xd[i]) > max) max = ABS(xd[i]);
643   }
644
645   return(max);
646 }
647
648 realtype N_VWrmsNorm_Serial(N_Vector x, N_Vector w)
649 {
650   long int i, N;
651   realtype sum, prodi, *xd, *wd;
652
653   sum = ZERO;
654   xd = wd = NULL;
655
656   N  = NV_LENGTH_S(x);
657   xd = NV_DATA_S(x);
658   wd = NV_DATA_S(w);
659
660   for (i = 0; i < N; i++) {
661     prodi = xd[i]*wd[i];
662     sum += SQR(prodi);
663   }
664
665   return(RSqrt(sum/N));
666 }
667
668 realtype N_VWrmsNormMask_Serial(N_Vector x, N_Vector w, N_Vector id)
669 {
670   long int i, N;
671   realtype sum, prodi, *xd, *wd, *idd;
672
673   sum = ZERO;
674   xd = wd = idd = NULL;
675
676   N  = NV_LENGTH_S(x);
677   xd  = NV_DATA_S(x);
678   wd  = NV_DATA_S(w);
679   idd = NV_DATA_S(id);
680
681   for (i = 0; i < N; i++) {
682     if (idd[i] > ZERO) {
683       prodi = xd[i]*wd[i];
684       sum += SQR(prodi);
685     }
686   }
687
688   return(RSqrt(sum / N));
689 }
690
691 realtype N_VMin_Serial(N_Vector x)
692 {
693   long int i, N;
694   realtype min, *xd;
695
696   xd = NULL;
697
698   N  = NV_LENGTH_S(x);
699   xd = NV_DATA_S(x);
700
701   min = xd[0];
702
703   for (i = 1; i < N; i++) {
704     if (xd[i] < min) min = xd[i];
705   }
706
707   return(min);
708 }
709
710 realtype N_VWL2Norm_Serial(N_Vector x, N_Vector w)
711 {
712   long int i, N;
713   realtype sum, prodi, *xd, *wd;
714
715   sum = ZERO;
716   xd = wd = NULL;
717
718   N  = NV_LENGTH_S(x);
719   xd = NV_DATA_S(x);
720   wd = NV_DATA_S(w);
721
722   for (i = 0; i < N; i++) {
723     prodi = xd[i]*wd[i];
724     sum += SQR(prodi);
725   }
726
727   return(RSqrt(sum));
728 }
729
730 realtype N_VL1Norm_Serial(N_Vector x)
731 {
732   long int i, N;
733   realtype sum, *xd;
734
735   sum = ZERO;
736   xd = NULL;
737
738   N  = NV_LENGTH_S(x);
739   xd = NV_DATA_S(x);
740   
741   for (i = 0; i<N; i++)  
742     sum += ABS(xd[i]);
743
744   return(sum);
745 }
746
747 void N_VCompare_Serial(realtype c, N_Vector x, N_Vector z)
748 {
749   long int i, N;
750   realtype *xd, *zd;
751
752   xd = zd = NULL;
753
754   N  = NV_LENGTH_S(x);
755   xd = NV_DATA_S(x);
756   zd = NV_DATA_S(z);
757
758   for (i = 0; i < N; i++) {
759     zd[i] = (ABS(xd[i]) >= c) ? ONE : ZERO;
760   }
761
762   return;
763 }
764
765 booleantype N_VInvTest_Serial(N_Vector x, N_Vector z)
766 {
767   long int i, N;
768   realtype *xd, *zd;
769
770   xd = zd = NULL;
771
772   N  = NV_LENGTH_S(x);
773   xd = NV_DATA_S(x);
774   zd = NV_DATA_S(z);
775
776   for (i = 0; i < N; i++) {
777     if (xd[i] == ZERO) return(FALSE);
778     zd[i] = ONE/xd[i];
779   }
780
781   return(TRUE);
782 }
783
784 booleantype N_VConstrMask_Serial(N_Vector c, N_Vector x, N_Vector m)
785 {
786   long int i, N;
787   booleantype test;
788   realtype *cd, *xd, *md;
789
790   cd = xd = md = NULL;
791
792   N  = NV_LENGTH_S(x);
793   xd = NV_DATA_S(x);
794   cd = NV_DATA_S(c);
795   md = NV_DATA_S(m);
796
797   test = TRUE;
798
799   for (i = 0; i < N; i++) {
800     md[i] = ZERO;
801     if (cd[i] == ZERO) continue;
802     if (cd[i] > ONEPT5 || cd[i] < -ONEPT5) {
803       if ( xd[i]*cd[i] <= ZERO) { test = FALSE; md[i] = ONE; }
804       continue;
805     }
806     if ( cd[i] > HALF || cd[i] < -HALF) {
807       if (xd[i]*cd[i] < ZERO ) { test = FALSE; md[i] = ONE; }
808     }
809   }
810
811   return(test);
812 }
813
814 realtype N_VMinQuotient_Serial(N_Vector num, N_Vector denom)
815 {
816   booleantype notEvenOnce;
817   long int i, N;
818   realtype *nd, *dd, min;
819
820   nd = dd = NULL;
821
822   N  = NV_LENGTH_S(num);
823   nd = NV_DATA_S(num);
824   dd = NV_DATA_S(denom);
825
826   notEvenOnce = TRUE;
827   min = BIG_REAL;
828
829   for (i = 0; i < N; i++) {
830     if (dd[i] == ZERO) continue;
831     else {
832       if (!notEvenOnce) min = MIN(min, nd[i]/dd[i]);
833       else {
834         min = nd[i]/dd[i];
835         notEvenOnce = FALSE;
836       }
837     }
838   }
839
840   return(min);
841 }
842
843 /*
844  * -----------------------------------------------------------------
845  * private functions
846  * -----------------------------------------------------------------
847  */
848
849 static void VCopy_Serial(N_Vector x, N_Vector z)
850 {
851   long int i, N;
852   realtype *xd, *zd;
853
854   xd = zd = NULL;
855
856   N  = NV_LENGTH_S(x);
857   xd = NV_DATA_S(x);
858   zd = NV_DATA_S(z);
859
860   for (i = 0; i < N; i++)
861     zd[i] = xd[i]; 
862
863   return;
864 }
865
866 static void VSum_Serial(N_Vector x, N_Vector y, N_Vector z)
867 {
868   long int i, N;
869   realtype *xd, *yd, *zd;
870
871   xd = yd = zd = NULL;
872
873   N  = NV_LENGTH_S(x);
874   xd = NV_DATA_S(x);
875   yd = NV_DATA_S(y);
876   zd = NV_DATA_S(z);
877
878   for (i = 0; i < N; i++)
879     zd[i] = xd[i]+yd[i];
880
881   return;
882 }
883
884 static void VDiff_Serial(N_Vector x, N_Vector y, N_Vector z)
885 {
886   long int i, N;
887   realtype *xd, *yd, *zd;
888
889   xd = yd = zd = NULL;
890
891   N  = NV_LENGTH_S(x);
892   xd = NV_DATA_S(x);
893   yd = NV_DATA_S(y);
894   zd = NV_DATA_S(z);
895
896   for (i = 0; i < N; i++)
897     zd[i] = xd[i]-yd[i];
898
899   return;
900 }
901
902 static void VNeg_Serial(N_Vector x, N_Vector z)
903 {
904   long int i, N;
905   realtype *xd, *zd;
906
907   xd = zd = NULL;
908
909   N  = NV_LENGTH_S(x);
910   xd = NV_DATA_S(x);
911   zd = NV_DATA_S(z);
912   
913   for (i = 0; i < N; i++)
914     zd[i] = -xd[i];
915
916   return;
917 }
918
919 static void VScaleSum_Serial(realtype c, N_Vector x, N_Vector y, N_Vector z)
920 {
921   long int i, N;
922   realtype *xd, *yd, *zd;
923
924   xd = yd = zd = NULL;
925
926   N  = NV_LENGTH_S(x);
927   xd = NV_DATA_S(x);
928   yd = NV_DATA_S(y);
929   zd = NV_DATA_S(z);
930
931   for (i = 0; i < N; i++)
932     zd[i] = c*(xd[i]+yd[i]);
933
934   return;
935 }
936
937 static void VScaleDiff_Serial(realtype c, N_Vector x, N_Vector y, N_Vector z)
938 {
939   long int i, N;
940   realtype *xd, *yd, *zd;
941
942   xd = yd = zd = NULL;
943
944   N  = NV_LENGTH_S(x);
945   xd = NV_DATA_S(x);
946   yd = NV_DATA_S(y);
947   zd = NV_DATA_S(z);
948
949   for (i = 0; i < N; i++)
950     zd[i] = c*(xd[i]-yd[i]);
951
952   return;
953 }
954
955 static void VLin1_Serial(realtype a, N_Vector x, N_Vector y, N_Vector z)
956 {
957   long int i, N;
958   realtype *xd, *yd, *zd;
959
960   xd = yd = zd = NULL;
961
962   N  = NV_LENGTH_S(x);
963   xd = NV_DATA_S(x);
964   yd = NV_DATA_S(y);
965   zd = NV_DATA_S(z);
966
967   for (i = 0; i < N; i++)
968     zd[i] = (a*xd[i])+yd[i];
969
970   return;
971 }
972
973 static void VLin2_Serial(realtype a, N_Vector x, N_Vector y, N_Vector z)
974 {
975   long int i, N;
976   realtype *xd, *yd, *zd;
977
978   xd = yd = zd = NULL;
979
980   N  = NV_LENGTH_S(x);
981   xd = NV_DATA_S(x);
982   yd = NV_DATA_S(y);
983   zd = NV_DATA_S(z);
984
985   for (i = 0; i < N; i++)
986     zd[i] = (a*xd[i])-yd[i];
987
988   return;
989 }
990
991 static void Vaxpy_Serial(realtype a, N_Vector x, N_Vector y)
992 {
993   long int i, N;
994   realtype *xd, *yd;
995
996   xd = yd = NULL;
997
998   N  = NV_LENGTH_S(x);
999   xd = NV_DATA_S(x);
1000   yd = NV_DATA_S(y);
1001
1002   if (a == ONE) {
1003     for (i = 0; i < N; i++)
1004       yd[i] += xd[i];
1005     return;
1006   }
1007
1008   if (a == -ONE) {
1009     for (i = 0; i < N; i++)
1010       yd[i] -= xd[i];
1011     return;
1012   }    
1013
1014   for (i = 0; i < N; i++)
1015     yd[i] += a*xd[i];
1016
1017   return;
1018 }
1019
1020 static void VScaleBy_Serial(realtype a, N_Vector x)
1021 {
1022   long int i, N;
1023   realtype *xd;
1024
1025   xd = NULL;
1026
1027   N  = NV_LENGTH_S(x);
1028   xd = NV_DATA_S(x);
1029
1030   for (i = 0; i < N; i++)
1031     xd[i] *= a;
1032
1033   return;
1034 }