/[ascend]/trunk/minos54/mi15blas.f
ViewVC logotype

Contents of /trunk/minos54/mi15blas.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (17 years, 9 months ago) by aw0a
File size: 18877 byte(s)
Setting up web subdirectory in repository
1 *********************************************************************
2 *
3 * File mi15blas fortran.
4 *
5 * dasum daxpy dcopy ddot dnrm2 dscal idamax
6 *
7 * These correspond to members of the BLAS package (Basic Linear
8 * Algebra Subprograms, Lawson et al. (1979), ACM TOMS 5, 3).
9 * If possible, they should be replaced by authentic BLAS routines
10 * tuned to the machine being used.
11 *
12 * The following utilities are also used:
13 *
14 * dddiv ddscl dload dnorm1
15 * hcopy hload icopy iload iload1
16 *
17 * These could be tuned to the machine being used.
18 * dload is used the most.
19 *
20 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21
22 *** dasum thru idamax taken
23 *** from netlib, Thu May 16 21:00:13 EDT 1991 ***
24 *** Declarations of the form dx(1) changed to dx(*)
25
26 double precision function dasum(n,dx,incx)
27 c
28 c takes the sum of the absolute values.
29 c jack dongarra, linpack, 3/11/78.
30 c
31 double precision dx(*),dtemp
32 integer i,incx,m,mp1,n,nincx
33 c
34 dasum = 0.0d0
35 dtemp = 0.0d0
36 if(n.le.0)return
37 if(incx.eq.1)go to 20
38 c
39 c code for increment not equal to 1
40 c
41 nincx = n*incx
42 do 10 i = 1,nincx,incx
43 dtemp = dtemp + dabs(dx(i))
44 10 continue
45 dasum = dtemp
46 return
47 c
48 c code for increment equal to 1
49 c
50 c
51 c clean-up loop
52 c
53 20 m = mod(n,6)
54 if( m .eq. 0 ) go to 40
55 do 30 i = 1,m
56 dtemp = dtemp + dabs(dx(i))
57 30 continue
58 if( n .lt. 6 ) go to 60
59 40 mp1 = m + 1
60 do 50 i = mp1,n,6
61 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2))
62 * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5))
63 50 continue
64 60 dasum = dtemp
65 return
66 end
67
68 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
69
70 subroutine daxpy(n,da,dx,incx,dy,incy)
71 c
72 c constant times a vector plus a vector.
73 c uses unrolled loops for increments equal to one.
74 c jack dongarra, linpack, 3/11/78.
75 c
76 double precision dx(*),dy(*),da
77 integer i,incx,incy,ix,iy,m,mp1,n
78 c
79 if(n.le.0)return
80 if (da .eq. 0.0d0) return
81 if(incx.eq.1.and.incy.eq.1)go to 20
82 c
83 c code for unequal increments or equal increments
84 c not equal to 1
85 c
86 ix = 1
87 iy = 1
88 if(incx.lt.0)ix = (-n+1)*incx + 1
89 if(incy.lt.0)iy = (-n+1)*incy + 1
90 do 10 i = 1,n
91 dy(iy) = dy(iy) + da*dx(ix)
92 ix = ix + incx
93 iy = iy + incy
94 10 continue
95 return
96 c
97 c code for both increments equal to 1
98 c
99 c
100 c clean-up loop
101 c
102 20 m = mod(n,4)
103 if( m .eq. 0 ) go to 40
104 do 30 i = 1,m
105 dy(i) = dy(i) + da*dx(i)
106 30 continue
107 if( n .lt. 4 ) return
108 40 mp1 = m + 1
109 do 50 i = mp1,n,4
110 dy(i) = dy(i) + da*dx(i)
111 dy(i + 1) = dy(i + 1) + da*dx(i + 1)
112 dy(i + 2) = dy(i + 2) + da*dx(i + 2)
113 dy(i + 3) = dy(i + 3) + da*dx(i + 3)
114 50 continue
115 return
116 end
117
118 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
119
120 subroutine dcopy(n,dx,incx,dy,incy)
121 c
122 c copies a vector, x, to a vector, y.
123 c uses unrolled loops for increments equal to one.
124 c jack dongarra, linpack, 3/11/78.
125 c
126 double precision dx(*),dy(*)
127 integer i,incx,incy,ix,iy,m,mp1,n
128 c
129 if(n.le.0)return
130 if(incx.eq.1.and.incy.eq.1)go to 20
131 c
132 c code for unequal increments or equal increments
133 c not equal to 1
134 c
135 ix = 1
136 iy = 1
137 if(incx.lt.0)ix = (-n+1)*incx + 1
138 if(incy.lt.0)iy = (-n+1)*incy + 1
139 do 10 i = 1,n
140 dy(iy) = dx(ix)
141 ix = ix + incx
142 iy = iy + incy
143 10 continue
144 return
145 c
146 c code for both increments equal to 1
147 c
148 c
149 c clean-up loop
150 c
151 20 m = mod(n,7)
152 if( m .eq. 0 ) go to 40
153 do 30 i = 1,m
154 dy(i) = dx(i)
155 30 continue
156 if( n .lt. 7 ) return
157 40 mp1 = m + 1
158 do 50 i = mp1,n,7
159 dy(i) = dx(i)
160 dy(i + 1) = dx(i + 1)
161 dy(i + 2) = dx(i + 2)
162 dy(i + 3) = dx(i + 3)
163 dy(i + 4) = dx(i + 4)
164 dy(i + 5) = dx(i + 5)
165 dy(i + 6) = dx(i + 6)
166 50 continue
167 return
168 end
169
170 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
171
172 double precision function ddot(n,dx,incx,dy,incy)
173 c
174 c forms the dot product of two vectors.
175 c uses unrolled loops for increments equal to one.
176 c jack dongarra, linpack, 3/11/78.
177 c
178 double precision dx(*),dy(*),dtemp
179 integer i,incx,incy,ix,iy,m,mp1,n
180 c
181 ddot = 0.0d0
182 dtemp = 0.0d0
183 if(n.le.0)return
184 if(incx.eq.1.and.incy.eq.1)go to 20
185 c
186 c code for unequal increments or equal increments
187 c not equal to 1
188 c
189 ix = 1
190 iy = 1
191 if(incx.lt.0)ix = (-n+1)*incx + 1
192 if(incy.lt.0)iy = (-n+1)*incy + 1
193 do 10 i = 1,n
194 dtemp = dtemp + dx(ix)*dy(iy)
195 ix = ix + incx
196 iy = iy + incy
197 10 continue
198 ddot = dtemp
199 return
200 c
201 c code for both increments equal to 1
202 c
203 c
204 c clean-up loop
205 c
206 20 m = mod(n,5)
207 if( m .eq. 0 ) go to 40
208 do 30 i = 1,m
209 dtemp = dtemp + dx(i)*dy(i)
210 30 continue
211 if( n .lt. 5 ) go to 60
212 40 mp1 = m + 1
213 do 50 i = mp1,n,5
214 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
215 * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
216 50 continue
217 60 ddot = dtemp
218 return
219 end
220
221 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
222
223 double precision function dnrm2 ( n, dx, incx)
224 integer next
225 double precision dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one
226 data zero, one /0.0d0, 1.0d0/
227 c
228 c euclidean norm of the n-vector stored in dx() with storage
229 c increment incx .
230 c if n .le. 0 return with result = 0.
231 c if n .ge. 1 then incx must be .ge. 1
232 c
233 c c.l.lawson, 1978 jan 08
234 c
235 c four phase method using two built-in constants that are
236 c hopefully applicable to all machines.
237 c cutlo = maximum of dsqrt(u/eps) over all known machines.
238 c cuthi = minimum of dsqrt(v) over all known machines.
239 c where
240 c eps = smallest no. such that eps + 1. .gt. 1.
241 c u = smallest positive no. (underflow limit)
242 c v = largest no. (overflow limit)
243 c
244 c brief outline of algorithm..
245 c
246 c phase 1 scans zero components.
247 c move to phase 2 when a component is nonzero and .le. cutlo
248 c move to phase 3 when a component is .gt. cutlo
249 c move to phase 4 when a component is .ge. cuthi/m
250 c where m = n for x() real and m = 2*n for complex.
251 c
252 c values for cutlo and cuthi..
253 c from the environmental parameters listed in the imsl converter
254 c document the limiting values are as follows..
255 c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
256 c univac and dec at 2**(-103)
257 c thus cutlo = 2**(-51) = 4.44089e-16
258 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
259 c thus cuthi = 2**(63.5) = 1.30438e19
260 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
261 c thus cutlo = 2**(-33.5) = 8.23181d-11
262 c cuthi, d.p. same as s.p. cuthi = 1.30438d19
263 c data cutlo, cuthi / 8.232d-11, 1.304d19 /
264 c data cutlo, cuthi / 4.441e-16, 1.304e19 /
265 data cutlo, cuthi / 8.232d-11, 1.304d19 /
266 c
267 if(n .gt. 0) go to 10
268 dnrm2 = zero
269 go to 300
270 c
271 10 assign 30 to next
272 sum = zero
273 nn = n * incx
274 c begin main loop
275 i = 1
276 20 go to next,(30, 50, 70, 110)
277 30 if( dabs(dx(i)) .gt. cutlo) go to 85
278 assign 50 to next
279 xmax = zero
280 c
281 c phase 1. sum is zero
282 c
283 50 if( dx(i) .eq. zero) go to 200
284 if( dabs(dx(i)) .gt. cutlo) go to 85
285 c
286 c prepare for phase 2.
287 assign 70 to next
288 go to 105
289 c
290 c prepare for phase 4.
291 c
292 100 i = j
293 assign 110 to next
294 sum = (sum / dx(i)) / dx(i)
295 105 xmax = dabs(dx(i))
296 go to 115
297 c
298 c phase 2. sum is small.
299 c scale to avoid destructive underflow.
300 c
301 70 if( dabs(dx(i)) .gt. cutlo ) go to 75
302 c
303 c common code for phases 2 and 4.
304 c in phase 4 sum is large. scale to avoid overflow.
305 c
306 110 if( dabs(dx(i)) .le. xmax ) go to 115
307 sum = one + sum * (xmax / dx(i))**2
308 xmax = dabs(dx(i))
309 go to 200
310 c
311 115 sum = sum + (dx(i)/xmax)**2
312 go to 200
313 c
314 c
315 c prepare for phase 3.
316 c
317 75 sum = (sum * xmax) * xmax
318 c
319 c
320 c for real or d.p. set hitest = cuthi/n
321 c for complex set hitest = cuthi/(2*n)
322 c
323 85 hitest = cuthi/float( n )
324 c
325 c phase 3. sum is mid-range. no scaling.
326 c
327 do 95 j =i,nn,incx
328 if(dabs(dx(j)) .ge. hitest) go to 100
329 95 sum = sum + dx(j)**2
330 dnrm2 = dsqrt( sum )
331 go to 300
332 c
333 200 continue
334 i = i + incx
335 if ( i .le. nn ) go to 20
336 c
337 c end of main loop.
338 c
339 c compute square root and adjust for scaling.
340 c
341 dnrm2 = xmax * dsqrt(sum)
342 300 continue
343 return
344 end
345
346 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
347
348 subroutine dscal(n,da,dx,incx)
349 c
350 c scales a vector by a constant.
351 c uses unrolled loops for increment equal to one.
352 c jack dongarra, linpack, 3/11/78.
353 c
354 double precision da,dx(*)
355 integer i,incx,m,mp1,n,nincx
356 c
357 if(n.le.0)return
358 if(incx.eq.1)go to 20
359 c
360 c code for increment not equal to 1
361 c
362 nincx = n*incx
363 do 10 i = 1,nincx,incx
364 dx(i) = da*dx(i)
365 10 continue
366 return
367 c
368 c code for increment equal to 1
369 c
370 c
371 c clean-up loop
372 c
373 20 m = mod(n,5)
374 if( m .eq. 0 ) go to 40
375 do 30 i = 1,m
376 dx(i) = da*dx(i)
377 30 continue
378 if( n .lt. 5 ) return
379 40 mp1 = m + 1
380 do 50 i = mp1,n,5
381 dx(i) = da*dx(i)
382 dx(i + 1) = da*dx(i + 1)
383 dx(i + 2) = da*dx(i + 2)
384 dx(i + 3) = da*dx(i + 3)
385 dx(i + 4) = da*dx(i + 4)
386 50 continue
387 return
388 end
389
390 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
391
392 integer function idamax(n,dx,incx)
393 c
394 c finds the index of element having max. absolute value.
395 c jack dongarra, linpack, 3/11/78.
396 c
397 double precision dx(*),dmax
398 integer i,incx,ix,n
399 c
400 idamax = 0
401 if( n .lt. 1 ) return
402 idamax = 1
403 if(n.eq.1)return
404 if(incx.eq.1)go to 20
405 c
406 c code for increment not equal to 1
407 c
408 ix = 1
409 dmax = dabs(dx(1))
410 ix = ix + incx
411 do 10 i = 2,n
412 if(dabs(dx(ix)).le.dmax) go to 5
413 idamax = i
414 dmax = dabs(dx(ix))
415 5 ix = ix + incx
416 10 continue
417 return
418 c
419 c code for increment equal to 1
420 c
421 20 dmax = dabs(dx(1))
422 do 30 i = 2,n
423 if(dabs(dx(i)).le.dmax) go to 30
424 idamax = i
425 dmax = dabs(dx(i))
426 30 continue
427 return
428 end
429
430 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
431
432 subroutine dddiv ( n, d, incd, x, incx )
433
434 implicit double precision (a-h,o-z)
435 double precision d(*), x(*)
436
437 * dddiv performs the diagonal scaling x = x / d.
438
439 integer i, id, ix
440 external dscal
441 intrinsic abs
442 parameter ( one = 1.0d+0 )
443
444 if (n .gt. 0) then
445 if (incd .eq. 0 .and. incx .ne. 0) then
446 call dscal ( n, one/d(1), x, abs(incx) )
447 else if (incd .eq. incx .and. incd .gt. 0) then
448 do 10 id = 1, 1 + (n - 1)*incd, incd
449 x(id) = x(id) / d(id)
450 10 continue
451 else
452 if (incx .ge. 0) then
453 ix = 1
454 else
455 ix = 1 - (n - 1)*incx
456 end if
457 if (incd .gt. 0) then
458 do 20 id = 1, 1 + (n - 1)*incd, incd
459 x(ix) = x(ix) / d(id)
460 ix = ix + incx
461 20 continue
462 else
463 id = 1 - (n - 1)*incd
464 do 30 i = 1, n
465 x(ix) = x(ix) / d(id)
466 id = id + incd
467 ix = ix + incx
468 30 continue
469 end if
470 end if
471 end if
472
473 * end of dddiv
474 end
475
476 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
477
478 subroutine ddscl ( n, d, incd, x, incx )
479
480 integer incd, incx, n
481 double precision d(*), x(*)
482
483 * ddscl performs the diagonal scaling x = d * x.
484
485 integer i, id, ix
486 external dscal
487 intrinsic abs
488
489 if (n .gt. 0) then
490 if (incd .eq. 0 .and. incx .ne. 0) then
491 call dscal ( n, d(1), x, abs(incx) )
492 else if (incd .eq. incx .and. incd .gt. 0) then
493 do 10 id = 1, 1 + (n - 1)*incd, incd
494 x(id) = d(id)*x(id)
495 10 continue
496 else
497 if (incx .ge. 0) then
498 ix = 1
499 else
500 ix = 1 - (n - 1)*incx
501 end if
502 if (incd .gt. 0) then
503 do 20 id = 1, 1 + (n - 1)*incd, incd
504 x(ix) = d(id)*x(ix)
505 ix = ix + incx
506 20 continue
507 else
508 id = 1 - (n - 1)*incd
509 do 30 i = 1, n
510 x(ix) = d(id)*x(ix)
511 id = id + incd
512 ix = ix + incx
513 30 continue
514 end if
515 end if
516 end if
517
518 * end of ddscl
519 end
520
521 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
522
523 subroutine dload ( n, const, x, incx )
524
525 double precision const
526 integer incx, n
527 double precision x(*)
528
529 * dload loads elements of x with const.
530
531 double precision zero
532 parameter ( zero = 0.0d+0 )
533 integer ix
534
535 if (n .gt. 0) then
536 if (incx .eq. 1 .and. const .eq. zero) then
537 do 10 ix = 1, n
538 x(ix) = zero
539 10 continue
540 else
541 do 20 ix = 1, 1 + (n - 1)*incx, incx
542 x(ix) = const
543 20 continue
544 end if
545 end if
546
547 * end of dload
548 end
549
550 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
551
552 function dnorm1( n, x, incx )
553
554 implicit double precision (a-h,o-z)
555 double precision x(*)
556
557 * dnorm1 returns the 1-norm of the vector x, scaled by root(n).
558 * This approximates an "average" element of x with some allowance
559 * for x being sparse.
560
561 intrinsic sqrt
562 external dasum
563
564 d = n
565 d = dasum ( n, x, incx ) / sqrt(d)
566 dnorm1 = d
567
568 * end of dnorm1
569 end
570
571 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
572
573 subroutine hcopy ( n, hx, incx, hy, incy )
574
575 integer*4 hx(*), hy(*)
576 integer incx, incy
577
578 * hcopy is the half-integer version of dcopy.
579 * In this version of MINOS we no longer use half integers.
580
581 integer ix, iy
582
583 if (n .gt. 0) then
584 if (incx .eq. incy .and. incy .gt. 0) then
585 do 10 iy = 1, 1 + (n - 1)*incy, incy
586 hy(iy) = hx(iy)
587 10 continue
588 else
589 if (incx .ge. 0) then
590 ix = 1
591 else
592 ix = 1 - (n - 1)*incx
593 end if
594 if (incy .gt. 0) then
595 do 20 iy = 1, 1 + ( n - 1 )*incy, incy
596 hy(iy) = hx(ix)
597 ix = ix + incx
598 20 continue
599 else
600 iy = 1 - (n - 1)*incy
601 do 30 i = 1, n
602 hy(iy) = hx(ix)
603 iy = iy + incy
604 ix = ix + incx
605 30 continue
606 end if
607 end if
608 end if
609
610 * end of hcopy
611 end
612
613 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
614
615 subroutine hload ( n, const, hx, incx )
616
617 integer incx, n
618 integer const
619 integer*4 hx(*)
620
621 * hload loads elements of hx with const.
622 * Beware that const is INTEGER, not half integer.
623
624 integer ix
625
626 if (n .gt. 0) then
627 if (incx .eq. 1 .and. const .eq. 0) then
628 do 10 ix = 1, n
629 hx(ix) = 0
630 10 continue
631 else
632 do 20 ix = 1, 1 + (n - 1)*incx, incx
633 hx(ix) = const
634 20 continue
635 end if
636 end if
637
638 * end of hload
639 end
640
641 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
642
643 subroutine icopy ( n, x, incx, y, incy )
644
645 integer x(*), y(*)
646 integer incx, incy
647
648 * icopy is the integer version of dcopy.
649
650 integer ix, iy
651
652 if (n .gt. 0) then
653 if (incx .eq. incy .and. incy .gt. 0) then
654 do 10 iy = 1, 1 + (n - 1)*incy, incy
655 y(iy) = x(iy)
656 10 continue
657 else
658 if (incx .ge. 0) then
659 ix = 1
660 else
661 ix = 1 - (n - 1)*incx
662 end if
663 if (incy .gt. 0) then
664 do 20 iy = 1, 1 + ( n - 1 )*incy, incy
665 y(iy) = x(ix)
666 ix = ix + incx
667 20 continue
668 else
669 iy = 1 - (n - 1)*incy
670 do 30 i = 1, n
671 y(iy) = x(ix)
672 iy = iy + incy
673 ix = ix + incx
674 30 continue
675 end if
676 end if
677 end if
678
679 * end of icopy
680 end
681
682 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
683
684 subroutine iload ( n, const, x, incx )
685
686 integer incx, n
687 integer const
688 integer x(*)
689
690 * iload loads elements of x with const.
691
692 integer ix
693
694 if (n .gt. 0) then
695 if (incx .eq. 1 .and. const .eq. 0) then
696 do 10 ix = 1, n
697 x(ix) = 0
698 10 continue
699 else
700 do 20 ix = 1, 1 + (n - 1)*incx, incx
701 x(ix) = const
702 20 continue
703 end if
704 end if
705
706 * end of iload
707 end
708
709 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
710
711 subroutine iload1( n, const, x, incx )
712
713 integer incx, n
714 integer const
715 integer x(*)
716
717 * iload1 loads elements of x with const, by calling iload.
718 * iload1 is needed in MINOS because iload is a file number.
719
720 call iload ( n, const, x, incx )
721
722 * end of iload1
723 end

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22