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 }