1
2 #include <stdlib.h>
3 #include <stdio.h>
4 /*#include <memory.h>*/
5 #include <string.h>
6 #include <math.h>
7 #include "myblas.h"
8
9 #ifdef FORTIFY
10 # include "lp_fortify.h"
11 #endif
12
13 /* ************************************************************************ */
14 /* Initialize BLAS interfacing routines */
15 /* ************************************************************************ */
16 MYBOOL mustinitBLAS = TRUE;
17 #if (defined WIN32) || (defined WIN64)
18 HINSTANCE hBLAS = NULL;
19 #else
20 void *hBLAS = NULL;
21 #endif
22
23
24 /* ************************************************************************ */
25 /* Function pointers for external BLAS library (C base 0) */
26 /* ************************************************************************ */
27 BLAS_dscal_func *BLAS_dscal;
28 BLAS_dcopy_func *BLAS_dcopy;
29 BLAS_daxpy_func *BLAS_daxpy;
30 BLAS_dswap_func *BLAS_dswap;
31 BLAS_ddot_func *BLAS_ddot;
32 BLAS_idamax_func *BLAS_idamax;
33 BLAS_idamin_func *BLAS_idamin;
34 BLAS_dload_func *BLAS_dload;
35 BLAS_dnormi_func *BLAS_dnormi;
36
37
38 /* ************************************************************************ */
39 /* Define the BLAS interfacing routines */
40 /* ************************************************************************ */
41
init_BLAS(void)42 void init_BLAS(void)
43 {
44 if(mustinitBLAS) {
45 load_BLAS(NULL);
46 mustinitBLAS = FALSE;
47 }
48 }
49
is_nativeBLAS(void)50 MYBOOL is_nativeBLAS(void)
51 {
52 #ifdef LoadableBlasLib
53 return( (MYBOOL) (hBLAS == NULL) );
54 #else
55 return( TRUE );
56 #endif
57 }
58
load_BLAS(char * libname)59 MYBOOL load_BLAS(char *libname)
60 {
61 MYBOOL result = TRUE;
62
63 #ifdef LoadableBlasLib
64 if(hBLAS != NULL) {
65 my_FreeLibrary(hBLAS);
66 }
67 #endif
68
69 if(libname == NULL) {
70 if(!mustinitBLAS && is_nativeBLAS())
71 return( FALSE );
72 BLAS_dscal = my_dscal;
73 BLAS_dcopy = my_dcopy;
74 BLAS_daxpy = my_daxpy;
75 BLAS_dswap = my_dswap;
76 BLAS_ddot = my_ddot;
77 BLAS_idamax = my_idamax;
78 BLAS_idamin = my_idamin;
79 BLAS_dload = my_dload;
80 BLAS_dnormi = my_dnormi;
81 if(mustinitBLAS)
82 mustinitBLAS = FALSE;
83 }
84 else {
85 #ifdef LoadableBlasLib
86 #if (defined WIN32) || (defined WIN64)
87 char *blasname = libname;
88 #else
89 /* First standardize UNIX .SO library name format. */
90 char blasname[260];
91 if(!so_stdname(blasname, libname, 260))
92 return( FALSE );
93 #endif
94 /* Get a handle to the Windows DLL module. */
95 hBLAS = my_LoadLibrary(blasname);
96
97 /* If the handle is valid, try to get the function addresses. */
98 result = (MYBOOL) (hBLAS != NULL);
99 if(result) {
100 BLAS_dscal = (BLAS_dscal_func *) my_GetProcAddress(hBLAS, BLAS_prec "scal");
101 BLAS_dcopy = (BLAS_dcopy_func *) my_GetProcAddress(hBLAS, BLAS_prec "copy");
102 BLAS_daxpy = (BLAS_daxpy_func *) my_GetProcAddress(hBLAS, BLAS_prec "axpy");
103 BLAS_dswap = (BLAS_dswap_func *) my_GetProcAddress(hBLAS, BLAS_prec "swap");
104 BLAS_ddot = (BLAS_ddot_func *) my_GetProcAddress(hBLAS, BLAS_prec "dot");
105 BLAS_idamax = (BLAS_idamax_func *) my_GetProcAddress(hBLAS, "i" BLAS_prec "amax");
106 BLAS_idamin = (BLAS_idamin_func *) my_GetProcAddress(hBLAS, "i" BLAS_prec "amin");
107 #if 0
108 BLAS_dload = (BLAS_dload_func *) my_GetProcAddress(hBLAS, BLAS_prec "load");
109 BLAS_dnormi = (BLAS_dnormi_func *) my_GetProcAddress(hBLAS, BLAS_prec "normi");
110 #endif
111 }
112 #endif
113 /* Do validation */
114 if(!result ||
115 ((BLAS_dscal == NULL) ||
116 (BLAS_dcopy == NULL) ||
117 (BLAS_daxpy == NULL) ||
118 (BLAS_dswap == NULL) ||
119 (BLAS_ddot == NULL) ||
120 (BLAS_idamax == NULL) ||
121 (BLAS_idamin == NULL) ||
122 (BLAS_dload == NULL) ||
123 (BLAS_dnormi == NULL))
124 ) {
125 load_BLAS(NULL);
126 result = FALSE;
127 }
128 }
129 return( result );
130 }
unload_BLAS(void)131 MYBOOL unload_BLAS(void)
132 {
133 return( load_BLAS(NULL) );
134 }
135
136
137 /* ************************************************************************ */
138 /* Now define the unoptimized local BLAS functions */
139 /* ************************************************************************ */
daxpylpsolve(int n,REAL da,REAL * dx,int incx,REAL * dy,int incy)140 void daxpylpsolve( int n, REAL da, REAL *dx, int incx, REAL *dy, int incy)
141 {
142 dx++;
143 dy++;
144 BLAS_daxpy( &n, &da, dx, &incx, dy, &incy);
145 }
my_daxpy(int * _n,REAL * _da,REAL * dx,int * _incx,REAL * dy,int * _incy)146 void BLAS_CALLMODEL my_daxpy( int *_n, REAL *_da, REAL *dx, int *_incx, REAL *dy, int *_incy)
147 {
148
149 /* constant times a vector plus a vector.
150 uses unrolled loops for increments equal to one.
151 jack dongarra, linpack, 3/11/78.
152 modified 12/3/93, array[1] declarations changed to array[*] */
153
154 int i, ix, iy;
155 #ifndef DOFASTMATH
156 int m, mp1;
157 #endif
158 register REAL rda;
159 REAL da = *_da;
160 int n = *_n, incx = *_incx, incy = *_incy;
161
162 if (n <= 0) return;
163 if (da == 0.0) return;
164
165 dx--;
166 dy--;
167 ix = 1;
168 iy = 1;
169 if (incx < 0)
170 ix = (-n+1)*incx + 1;
171 if (incy < 0)
172 iy = (-n+1)*incy + 1;
173 rda = da;
174
175 /* CPU intensive loop; option to do pointer arithmetic */
176 #if defined DOFASTMATH
177 {
178 REAL *xptr, *yptr;
179 for (i = 1, xptr = dx + ix, yptr = dy + iy;
180 i <= n; i++, xptr += incx, yptr += incy)
181 (*yptr) += rda*(*xptr);
182 }
183 #else
184
185 if (incx==1 && incy==1) goto x20;
186
187 /* code for unequal increments or equal increments not equal to 1 */
188 for (i = 1; i<=n; i++) {
189 dy[iy]+= rda*dx[ix];
190 ix+= incx;
191 iy+= incy;
192 }
193 return;
194
195 /* code for both increments equal to 1 */
196
197 /* clean-up loop */
198 x20:
199 m = n % 4;
200 if (m == 0) goto x40;
201 for (i = 1; i<=m; i++)
202 dy[i]+= rda*dx[i];
203 if(n < 4) return;
204 x40:
205 mp1 = m + 1;
206 for (i = mp1; i<=n; i=i+4) {
207 dy[i]+= rda*dx[i];
208 dy[i + 1]+= rda*dx[i + 1];
209 dy[i + 2]+= rda*dx[i + 2];
210 dy[i + 3]+= rda*dx[i + 3];
211 }
212 #endif
213 }
214
215
216 /* ************************************************************************ */
dcopylpsolve(int n,REAL * dx,int incx,REAL * dy,int incy)217 void dcopylpsolve( int n, REAL *dx, int incx, REAL *dy, int incy)
218 {
219 dx++;
220 dy++;
221 BLAS_dcopy( &n, dx, &incx, dy, &incy);
222 }
223
my_dcopy(int * _n,REAL * dx,int * _incx,REAL * dy,int * _incy)224 void BLAS_CALLMODEL my_dcopy (int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy)
225 {
226
227 /* copies a vector, x, to a vector, y.
228 uses unrolled loops for increments equal to one.
229 jack dongarra, linpack, 3/11/78.
230 modified 12/3/93, array[1] declarations changed to array[*] */
231
232 int i, ix, iy;
233 #ifndef DOFASTMATH
234 int m, mp1;
235 #endif
236 int n = *_n, incx = *_incx, incy = *_incy;
237
238 if(n<=0)
239 return;
240
241 dx--;
242 dy--;
243 ix = 1;
244 iy = 1;
245 if(incx<0)
246 ix = (-n+1)*incx + 1;
247 if(incy<0)
248 iy = (-n+1)*incy + 1;
249
250
251 /* CPU intensive loop; option to do pointer arithmetic */
252 #if defined DOFASTMATH
253 {
254 REAL *xptr, *yptr;
255 for (i = 1, xptr = dx + ix, yptr = dy + iy;
256 i <= n; i++, xptr += incx, yptr += incy)
257 (*yptr) = (*xptr);
258 }
259 #else
260
261 if(incx==1 && incy==1)
262 goto x20;
263
264 /* code for unequal increments or equal increments not equal to 1 */
265
266 for(i = 1; i<=n; i++) {
267 dy[iy] = dx[ix];
268 ix+= incx;
269 iy+= incy;
270 }
271 return;
272
273 /* code for both increments equal to 1 */
274
275 /* version with fast machine copy logic (requires memory.h or string.h) */
276 x20:
277 m = n % 7;
278 if (m == 0) goto x40;
279 for (i = 1; i<=m; i++)
280 dy[i] = dx[i];
281 if (n < 7) return;
282
283 x40:
284 mp1 = m + 1;
285 for (i = mp1; i<=n; i=i+7) {
286 dy[i] = dx[i];
287 dy[i + 1] = dx[i + 1];
288 dy[i + 2] = dx[i + 2];
289 dy[i + 3] = dx[i + 3];
290 dy[i + 4] = dx[i + 4];
291 dy[i + 5] = dx[i + 5];
292 dy[i + 6] = dx[i + 6];
293 }
294 #endif
295 }
296
297
298 /* ************************************************************************ */
299
dscallpsolve(int n,REAL da,REAL * dx,int incx)300 void dscallpsolve (int n, REAL da, REAL *dx, int incx)
301 {
302 dx++;
303 BLAS_dscal (&n, &da, dx, &incx);
304 }
305
my_dscal(int * _n,REAL * _da,REAL * dx,int * _incx)306 void BLAS_CALLMODEL my_dscal (int *_n, REAL *_da, REAL *dx, int *_incx)
307 {
308
309 /* Multiply a vector by a constant.
310
311 --Input--
312 N number of elements in input vector(s)
313 DA double precision scale factor
314 DX double precision vector with N elements
315 INCX storage spacing between elements of DX
316
317 --Output--
318 DX double precision result (unchanged if N.LE.0)
319
320 Replace double precision DX by double precision DA*DX.
321 For I = 0 to N-1, replace DX(IX+I*INCX) with DA * DX(IX+I*INCX),
322 where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX. */
323
324 int i;
325 #ifndef DOFASTMATH
326 int m, mp1, ix;
327 #endif
328 register REAL rda;
329 REAL da = *_da;
330 int n = *_n, incx = *_incx;
331
332 if (n <= 0)
333 return;
334 rda = da;
335
336 dx--;
337
338 /* Optionally do fast pointer arithmetic */
339 #if defined DOFASTMATH
340 {
341 REAL *xptr;
342 for (i = 1, xptr = dx + 1; i <= n; i++, xptr += incx)
343 (*xptr) *= rda;
344 }
345 #else
346
347 if (incx == 1)
348 goto x20;
349
350 /* Code for increment not equal to 1 */
351 ix = 1;
352 if (incx < 0)
353 ix = (-n+1)*incx + 1;
354 for(i = 1; i <= n; i++, ix += incx)
355 dx[ix] *= rda;
356 return;
357
358 /* Code for increment equal to 1. */
359 /* Clean-up loop so remaining vector length is a multiple of 5. */
360 x20:
361 m = n % 5;
362 if (m == 0) goto x40;
363 for( i = 1; i <= m; i++)
364 dx[i] *= rda;
365 if (n < 5)
366 return;
367 x40:
368 mp1 = m + 1;
369 for(i = mp1; i <= n; i += 5) {
370 dx[i] *= rda;
371 dx[i+1] *= rda;
372 dx[i+2] *= rda;
373 dx[i+3] *= rda;
374 dx[i+4] *= rda;
375 }
376 #endif
377 }
378
379
380 /* ************************************************************************ */
381
ddotlpsolve(int n,REAL * dx,int incx,REAL * dy,int incy)382 REAL ddotlpsolve(int n, REAL *dx, int incx, REAL *dy, int incy)
383 {
384 dx++;
385 dy++;
386 return( BLAS_ddot (&n, dx, &incx, dy, &incy) );
387 }
388
my_ddot(int * _n,REAL * dx,int * _incx,REAL * dy,int * _incy)389 REAL BLAS_CALLMODEL my_ddot(int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy)
390 {
391
392 /* forms the dot product of two vectors.
393 uses unrolled loops for increments equal to one.
394 jack dongarra, linpack, 3/11/78.
395 modified 12/3/93, array[1] declarations changed to array[*] */
396
397 register REAL dtemp;
398 int i, ix, iy;
399 #ifndef DOFASTMATH
400 int m, mp1;
401 #endif
402 int n = *_n, incx = *_incx, incy = *_incy;
403
404 dtemp = 0.0;
405 if (n<=0)
406 return( (REAL) dtemp);
407
408 dx--;
409 dy--;
410 ix = 1;
411 iy = 1;
412 if (incx<0)
413 ix = (-n+1)*incx + 1;
414 if (incy<0)
415 iy = (-n+1)*incy + 1;
416
417 /* CPU intensive loop; option to do pointer arithmetic */
418
419 #if defined DOFASTMATH
420 {
421 REAL *xptr, *yptr;
422 for (i = 1, xptr = dx + ix, yptr = dy + iy;
423 i <= n; i++, xptr += incx, yptr += incy)
424 dtemp+= (*yptr)*(*xptr);
425 }
426 #else
427
428 if (incx==1 && incy==1) goto x20;
429
430 /* code for unequal increments or equal increments not equal to 1 */
431
432 for (i = 1; i<=n; i++) {
433 dtemp+= dx[ix]*dy[iy];
434 ix+= incx;
435 iy+= incy;
436 }
437 return(dtemp);
438
439 /* code for both increments equal to 1 */
440
441 /* clean-up loop */
442
443 x20:
444 m = n % 5;
445 if (m == 0) goto x40;
446 for (i = 1; i<=m; i++)
447 dtemp+= dx[i]*dy[i];
448 if (n < 5) goto x60;
449
450 x40:
451 mp1 = m + 1;
452 for (i = mp1; i<=n; i=i+5)
453 dtemp+= dx[i]*dy[i] + dx[i + 1]*dy[i + 1] +
454 dx[i + 2]*dy[i + 2] + dx[i + 3]*dy[i + 3] + dx[i + 4]*dy[i + 4];
455
456 x60:
457 #endif
458 return(dtemp);
459 }
460
461
462 /* ************************************************************************ */
463
dswaplpsolve(int n,REAL * dx,int incx,REAL * dy,int incy)464 void dswaplpsolve( int n, REAL *dx, int incx, REAL *dy, int incy )
465 {
466 dx++;
467 dy++;
468 BLAS_dswap( &n, dx, &incx, dy, &incy );
469 }
470
my_dswap(int * _n,REAL * dx,int * _incx,REAL * dy,int * _incy)471 void BLAS_CALLMODEL my_dswap( int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy )
472 {
473 int i, ix, iy;
474 #ifndef DOFASTMATH
475 int m, mp1, ns;
476 REAL dtemp2, dtemp3;
477 #endif
478 register REAL dtemp1;
479 int n = *_n, incx = *_incx, incy = *_incy;
480
481 if (n <= 0) return;
482
483 dx--;
484 dy--;
485 ix = 1;
486 iy = 1;
487 if (incx < 0)
488 ix = (-n+1)*incx + 1;
489 if (incy < 0)
490 iy = (-n+1)*incy + 1;
491
492 /* CPU intensive loop; option to do pointer arithmetic */
493 #if defined DOFASTMATH
494 {
495 REAL *xptr, *yptr;
496 for (i = 1, xptr = dx + ix, yptr = dy + iy;
497 i <= n; i++, xptr += incx, yptr += incy) {
498 dtemp1 = (*xptr);
499 (*xptr) = (*yptr);
500 (*yptr) = dtemp1;
501 }
502 }
503 #else
504
505 if (incx == incy) {
506 if (incx <= 0) goto x5;
507 if (incx == 1) goto x20;
508 goto x60;
509 }
510
511 /* code for unequal or nonpositive increments. */
512 x5:
513 for (i = 1; i<=n; i++) {
514 dtemp1 = dx[ix];
515 dx[ix] = dy[iy];
516 dy[iy] = dtemp1;
517 ix+= incx;
518 iy+= incy;
519 }
520 return;
521
522 /* code for both increments equal to 1.
523 clean-up loop so remaining vector length is a multiple of 3. */
524 x20:
525 m = n % 3;
526 if (m == 0) goto x40;
527 for (i = 1; i<=m; i++) {
528 dtemp1 = dx[i];
529 dx[i] = dy[i];
530 dy[i] = dtemp1;
531 }
532 if (n < 3) return;
533
534 x40:
535 mp1 = m + 1;
536 for (i = mp1; i<=n; i=i+3) {
537 dtemp1 = dx[i];
538 dtemp2 = dx[i+1];
539 dtemp3 = dx[i+2];
540 dx[i] = dy[i];
541 dx[i+1] = dy[i+1];
542 dx[i+2] = dy[i+2];
543 dy[i] = dtemp1;
544 dy[i+1] = dtemp2;
545 dy[i+2] = dtemp3;
546 }
547 return;
548
549 /* code for equal, positive, non-unit increments. */
550 x60:
551 ns = n*incx;
552 for (i = 1; i<=ns; i=i+incx) {
553 dtemp1 = dx[i];
554 dx[i] = dy[i];
555 dy[i] = dtemp1;
556 }
557 #endif
558
559 }
560
561
562 /* ************************************************************************ */
563
dload(int n,REAL da,REAL * dx,int incx)564 void dload(int n, REAL da, REAL *dx, int incx)
565 {
566 dx++;
567 BLAS_dload (&n, &da, dx, &incx);
568 }
569
my_dload(int * _n,REAL * _da,REAL * dx,int * _incx)570 void BLAS_CALLMODEL my_dload (int *_n, REAL *_da, REAL *dx, int *_incx)
571 {
572 /* copies a scalar, a, to a vector, x.
573 uses unrolled loops when incx equals one.
574
575 To change the precision of this program, run the change
576 program on dload.f
577 Alternatively, to make a single precision version append a
578 comment character to the start of all lines between sequential
579 precision > double
580 and
581 end precision > double
582 comments and delete the comment character at the start of all
583 lines between sequential
584 precision > single
585 and
586 end precision > single
587 comments. To make a double precision version interchange
588 the append and delete operations in these instructions. */
589
590 int i, ix, m, mp1;
591 REAL da = *_da;
592 int n = *_n, incx = *_incx;
593
594 if (n<=0) return;
595 dx--;
596 if (incx==1) goto x20;
597
598 /* code for incx not equal to 1 */
599
600 ix = 1;
601 if (incx<0)
602 ix = (-n+1)*incx + 1;
603 for (i = 1; i<=n; i++) {
604 dx[ix] = da;
605 ix+= incx;
606 }
607 return;
608
609 /* code for incx equal to 1 and clean-up loop */
610
611 x20:
612 m = n % 7;
613 if (m == 0) goto x40;
614 for (i = 1; i<=m; i++)
615 dx[i] = da;
616 if (n < 7) return;
617
618 x40:
619 mp1 = m + 1;
620 for (i = mp1; i<=n; i=i+7) {
621 dx[i] = da;
622 dx[i + 1] = da;
623 dx[i + 2] = da;
624 dx[i + 3] = da;
625 dx[i + 4] = da;
626 dx[i + 5] = da;
627 dx[i + 6] = da;
628 }
629 }
630
631 /* ************************************************************************ */
idamaxlpsolve(int n,REAL * x,int is)632 int idamaxlpsolve( int n, REAL *x, int is )
633 {
634 x++;
635 return ( BLAS_idamax( &n, x, &is ) );
636 }
637
idaminlpsolve(int n,REAL * x,int is)638 int idaminlpsolve( int n, REAL *x, int is )
639 {
640 x++;
641 return ( BLAS_idamin( &n, x, &is ) );
642 }
643
my_idamax(int * _n,REAL * x,int * _is)644 int BLAS_CALLMODEL my_idamax( int *_n, REAL *x, int *_is )
645 {
646 register REAL xmax, xtest;
647 int i, imax = 0;
648 #if !defined DOFASTMATH
649 int ii;
650 #endif
651 int n = *_n, is = *_is;
652
653 if((n < 1) || (is <= 0))
654 return(imax);
655 imax = 1;
656 if(n == 1)
657 return(imax);
658
659 #if defined DOFASTMATH
660 xmax = fabs(*x);
661 for (i = 2, x += is; i <= n; i++, x += is) {
662 xtest = fabs(*x);
663 if(xtest > xmax) {
664 xmax = xtest;
665 imax = i;
666 }
667 }
668 #else
669 x--;
670 ii = 1;
671 xmax = fabs(x[ii]);
672 for(i = 2, ii+ = is; i <= n; i++, ii+ = is) {
673 xtest = fabs(x[ii]);
674 if(xtest > xmax) {
675 xmax = xtest;
676 imax = i;
677 }
678 }
679 #endif
680 return(imax);
681 }
682
my_idamin(int * _n,REAL * x,int * _is)683 int BLAS_CALLMODEL my_idamin( int *_n, REAL *x, int *_is )
684 {
685 register REAL xmin, xtest;
686 int i, imin = 0;
687 #if !defined DOFASTMATH
688 int ii;
689 #endif
690 int n = *_n, is = *_is;
691
692 if((n < 1) || (is <= 0))
693 return(imin);
694 imin = 1;
695 if(n == 1)
696 return(imin);
697
698 #if defined DOFASTMATH
699 xmin = fabs(*x);
700 for (i = 2, x += is; i <= n; i++, x += is) {
701 xtest = fabs(*x);
702 if(xtest < xmin) {
703 xmin = xtest;
704 imin = i;
705 }
706 }
707 #else
708 x--;
709 ii = 1;
710 xmin = fabs(x[ii]);
711 for(i = 2, ii+ = is; i <= n; i++, ii+ = is) {
712 xtest = fabs(x[ii]);
713 if(xtest < xmin) {
714 xmin = xtest;
715 imin = i;
716 }
717 }
718 #endif
719 return(imin);
720 }
721
722 /* ************************************************************************ */
dnormi(int n,REAL * x)723 REAL dnormi( int n, REAL *x )
724 {
725 x++;
726 return( BLAS_dnormi( &n, x ) );
727 }
728
my_dnormi(int * _n,REAL * x)729 REAL BLAS_CALLMODEL my_dnormi( int *_n, REAL *x )
730 {
731 /* ===============================================================
732 dnormi returns the infinity-norm of the vector x.
733 =============================================================== */
734 int j;
735 register REAL hold, absval;
736 int n = *_n;
737
738 x--;
739 hold = 0.0;
740 /* for(j = 1; j <= n; j++) */
741 for(j = n; j > 0; j--) {
742 absval = fabs(x[j]);
743 hold = MAX( hold, absval );
744 }
745
746 return( hold );
747 }
748
749
750 /* ************************************************************************ */
751 /* Subvector and submatrix access routines (Fortran compatibility) */
752 /* ************************************************************************ */
753
754 #ifndef UseMacroVector
subvec(int item)755 int subvec( int item)
756 {
757 return( item-1 );
758 }
759 #endif
760
submat(int nrowb,int row,int col)761 int submat( int nrowb, int row, int col)
762 {
763 return( nrowb*(col-1) + subvec(row) );
764 }
765
posmat(int nrowb,int row,int col)766 int posmat( int nrowb, int row, int col)
767 {
768 return( submat(nrowb, row, col)+BLAS_BASE );
769 }
770
771 /* ************************************************************************ */
772 /* Randomization functions */
773 /* ************************************************************************ */
774
randomseed(int seeds[])775 void randomseed(int seeds[])
776 /* Simply create some default seed values */
777 {
778 seeds[1] = 123456;
779 seeds[2] = 234567;
780 seeds[3] = 345678;
781 }
782
randomdens(int n,REAL * x,REAL r1,REAL r2,REAL densty,int * seeds)783 void randomdens( int n, REAL *x, REAL r1, REAL r2, REAL densty, int *seeds )
784 {
785 /* ------------------------------------------------------------------
786 random generates a vector x[*] of random numbers
787 in the range (r1, r2) with (approximate) specified density.
788 seeds[*] must be initialized before the first call.
789 ------------------------------------------------------------------ */
790
791 int i;
792 REAL *y;
793
794 y = (REAL *) malloc(sizeof(*y) * (n+1));
795 ddrand( n, x, 1, seeds );
796 ddrand( n, y, 1, seeds );
797
798 for (i = 1; i<=n; i++) {
799 if (y[i] < densty)
800 x[i] = r1 + (r2 - r1) * x[i];
801 else
802 x[i] = 0.0;
803 }
804 free(y);
805 }
806
807
808 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
809
ddrand(int n,REAL * x,int incx,int * seeds)810 void ddrand( int n, REAL *x, int incx, int *seeds )
811 {
812
813 /* ------------------------------------------------------------------
814 ddrand fills a vector x with uniformly distributed random numbers
815 in the interval (0, 1) using a method due to Wichman and Hill.
816
817 seeds[1..3] should be set to integer values
818 between 1 and 30000 before the first entry.
819
820 Integer arithmetic up to 30323 is required.
821
822 Blatantly copied from Wichman and Hill 19-January-1987.
823 14-Feb-94. Original version.
824 30 Jun 1999. seeds stored in an array.
825 30 Jun 1999. This version of ddrand.
826 ------------------------------------------------------------------ */
827
828 int ix, xix;
829
830 if (n < 1) return;
831
832 for (ix = 1; ix<=1+(n-1)*incx; ix=ix+incx) {
833 seeds[1] = 171*(seeds[1] % 177) - 2*(seeds[1]/177);
834 seeds[2] = 172*(seeds[2] % 176) - 35*(seeds[2]/176);
835 seeds[3] = 170*(seeds[3] % 178) - 63*(seeds[3]/178);
836
837 if (seeds[1] < 0) seeds[1] = seeds[1] + 30269;
838 if (seeds[2] < 0) seeds[2] = seeds[2] + 30307;
839 if (seeds[3] < 0) seeds[3] = seeds[3] + 30323;
840
841 x[ix] = ((REAL) seeds[1])/30269.0 +
842 ((REAL) seeds[2])/30307.0 +
843 ((REAL) seeds[3])/30323.0;
844 xix = (int) x[ix];
845 x[ix] = fabs(x[ix] - xix);
846 }
847
848 }
849
850