New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

lib_isml.f90 in trunk/NEMO/OPA_SRC – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/lib_isml.f90 @ 88

Last change on this file since 88 was 88, checked in by opalod, 20 years ago

CT : UPDATE057 : # General syntax, alignement, comments corrections

# l_ctl alone replace the set (l_ctl .AND. lwp)
# Add of diagnostics which are activated when using l_ctl logical

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 19.1 KB

Line
1	!
2	! subroutines for PCG or SOR solvers
3	! (used if the ISML library is not available)
4	!
5	! linrg
6	! gauss
7	! vmov
8	! desremopt
9	! dtrsv
10	! dger
11	! xerbla
12	! lsame
13	! folr (empty)
14	!
15	!---------------------------------------------------------
16	SUBROUTINE linrg(kn,pa,klda,painv,kldainv)
17
18	!! compute inverse matrix
19
20	IMPLICIT NONE
21	INTEGER kn,klda,kldainv
22	REAL (kind=8) :: pa(kn,kn),painv(kn,kn)
23	REAL (kind=8) :: zb(kn,kn)
24	REAL (kind=8) :: zv(kn)
25	INTEGER iplin(kn)
26	INTEGER ji
27
28	IF( kn /= klda .OR. kn /= kldainv ) THEN
29	write(0,*)'change your parameters'
30	STOP
31	ENDIF
32
33	CALL vmov( kn*kn, pa, painv )
34
35	CALL gauss( kn, painv, iplin, zv )
36
37	zb(:,:) = 0.e0
38	DO ji = 1, kn
39	zb(ji,ji) = 1.e0
40	CALL desremopt( kn, painv, iplin, zb(1,ji), zb(1,ji), zv )
41	END DO
42	CALL vmov( kn*kn, zb, painv )
43
44	END SUBROUTINE linrg
45	!---------------------------------------------------------
46	SUBROUTINE gauss(kn,pa,kplin,pv)
47
48	IMPLICIT NONE
49	INTEGER kn
50	INTEGER ji,jj,jk
51	INTEGER ik,ipp
52	REAL (kind=8) :: pa(kn,kn),pv(kn)
53	!! REAL (kind=8) :: zpivmax,zalpha
54	REAL (kind=8) :: zalpha
55	INTEGER kplin(kn)
56	INTEGER isamax
57	EXTERNAL isamax
58
59	! factorisation de Gauss de la matrice a avec pivot partiel .
60	! initialisation des pointeurs .
61	DO ji=1,kn
62	kplin(ji)=ji
63	END DO
64	DO jk=1,kn-1
65	! recherche du pivot maximal .
66	!! ik=jk
67	!! zpivmax=dabs(pa(jk,jk))
68	!! DO ji=jk,kn
69	!! IF(dabs(pa(ji,jk)) > zpivmax) THEN
70	!! zpivmax=dabs(pa(ji,jk))
71	!! ik=ji
72	!! ENDIF
73	!! END DO
74	ik=isamax( kn-jk+1, pa(jk,jk) )+jk-1
75	! permutation de la ligne jk et de la ligne ik .
76	IF(jk == 58) THEN
77	PRINT *,'matrix ',(pa(jk,ji),ji=1,kn)
78	PRINT *,' pivot ',ik,kplin(ik),kplin(jk)
79	ENDIF
80	ipp=kplin(ik)
81	kplin(ik)=kplin(jk)
82	kplin(jk)=ipp
83	DO jj=1,kn
84	pv(jj)=pa(ik,jj)
85	pa(ik,jj)=pa(jk,jj)
86	pa(jk,jj)=pv(jj)
87	END DO
88	IF(jk == 58) THEN
89	PRINT *,'matrix ',(pa(jk,ji),ji=1,kn)
90	PRINT *,' pivot ',ik,kplin(ik),kplin(jk)
91	ENDIF
92	! calcul des coefficients de la colonne k ligne a ligne .
93	DO ji=jk+1,kn
94	IF(pa(jk,jk) == 0) THEN
95	PRINT *,'probleme diagonale nulle',jk,pa(jk,jk)
96	pa(ji,jk)=pa(ji,jk)/1.E-20
97	ENDIF
98	IF(pa(jk,jk) /= 0) pa(ji,jk)=pa(ji,jk)/pa(jk,jk)
99	END DO
100	!! DO ji=jk+1,kn
101	!! DO jj=jk+1,kn
102	!! pa(ji,jj)=pa(ji,jj)-pa(ji,jk)*pa(jk,jj)
103	!! END DO
104	!! END DO
105	zalpha=-1.
106	CALL dger(kn-jk,kn-jk,zalpha,pa(jk+1,jk),1,pa(jk,jk+1),kn, &
107	pa(jk+1,jk+1),kn)
108	END DO
109
110	END SUBROUTINE gauss
111	!---------------------------------------------------------
112	FUNCTION isamax( I, X )
113	DIMENSION X(I)
114	ISAMAX = 0
115	XMIN = -1e+50
116	DO N = 1, I
117	IF(ABS(X(N)) > XMIN ) THEN
118	XMIN = X(N)
119	ISAMAX = N
120	ENDIF
121	END DO
122	END FUNCTION isamax
123	!---------------------------------------------------------
124	SUBROUTINE vmov(kn,px,py)
125
126	IMPLICIT NONE
127	INTEGER kn
128	REAL (kind=8) :: px(kn),py(kn)
129	INTEGER ji
130
131	DO ji=1,kn
132	py(ji)=px(ji)
133	END DO
134
135	END SUBROUTINE vmov
136	!---------------------------------------------------------
137	subroutine desremopt(n,a,plin,y,x,v)
138	implicit none
139	integer n,i, j0
140	!! integer n,i,j,j0
141	real (kind=8) :: a(n,n),x(n),y(n),v(n)
142	integer plin(n)
143	! descente remontee du systeme .
144	! initialisation du vecteur resultat .
145	! prise en compte de la permutation des lignes .
146	do i=1,n
147	v(i)=y(plin(i))
148	end do
149	do i=1,n
150	if(v(i) /= 0.) then
151	j0=i-1
152	goto 1
153	endif
154	end do
155	1 continue
156	! descente du systeme L v = v , L est a diagonale unitaire .
157	!! do j=j0+1,n
158	!! do i=j+1,n
159	!! v(i)=v(i)-a(i,j)*v(j)
160	!! end do
161	!! end do
162	call dtrsv('L','N','U',n-j0,a(j0+1,j0+1),n,v(j0+1),1)
163	! remontee du systeme U v = v .
164	!! do j=n,1,-1
165	!! v(j)=v(j)/a(j,j)
166	!! do i=1,j-1
167	!! v(i)=v(i)-a(i,j)*v(j)
168	!! end do
169	!! end do
170	call dtrsv('U','N','N',n,a,n,v,1)
171	! prise en compte de la permutation des colonnes .
172	do i=1,n
173	x(i)=v(i)
174	end do
175
176	end SUBROUTINE desremopt
177	!---------------------------------------------------------
178	SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
179	!! .. Scalar Arguments ..
180	INTEGER INCX, LDA, N
181	CHARACTER (len=1) :: DIAG, TRANS, UPLO
182	!! .. Array Arguments ..
183	! DOUBLE PRECISION A( LDA, * ), X( * )
184	REAL (kind=8) :: A( LDA, * ), X( * )
185	!! ..
186
187	!! Purpose
188	!! =======
189
190	!! DTRSV solves one of the systems of equations
191
192	!! Ax = b, or A'x = b,
193
194	!! where b and x are n element vectors and A is an n by n unit, or
195	!! non-unit, upper or lower triangular matrix.
196
197	!! No test for singularity or near-singularity is included in this
198	!! routine. Such tests must be performed before calling this routine.
199
200	!! Parameters
201	!! ==========
202
203	!! UPLO - CHARACTER*1.
204	!! On entry, UPLO specifies whether the matrix is an upper or
205	!! lower triangular matrix as follows:
206
207	!! UPLO = 'U' or 'u' A is an upper triangular matrix.
208
209	!! UPLO = 'L' or 'l' A is a lower triangular matrix.
210
211	!! Unchanged on exit.
212
213	!! TRANS - CHARACTER*1.
214	!! On entry, TRANS specifies the equations to be solved as
215	!! follows:
216
217	!! TRANS = 'N' or 'n' A*x = b.
218
219	!! TRANS = 'T' or 't' A'*x = b.
220
221	!! TRANS = 'C' or 'c' A'*x = b.
222
223	!! Unchanged on exit.
224
225	!! DIAG - CHARACTER*1.
226	!! On entry, DIAG specifies whether or not A is unit
227	!! triangular as follows:
228
229	!! DIAG = 'U' or 'u' A is assumed to be unit triangular.
230
231	!! DIAG = 'N' or 'n' A is not assumed to be unit
232	!! triangular.
233
234	!! Unchanged on exit.
235
236	!! N - INTEGER.
237	!! On entry, N specifies the order of the matrix A.
238	!! N must be at least zero.
239	!! Unchanged on exit.
240
241	!! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
242	!! Before entry with UPLO = 'U' or 'u', the leading n by n
243	!! upper triangular part of the array A must contain the upper
244	!! triangular matrix and the strictly lower triangular part of
245	!! A is not referenced.
246	!! Before entry with UPLO = 'L' or 'l', the leading n by n
247	!! lower triangular part of the array A must contain the lower
248	!! triangular matrix and the strictly upper triangular part of
249	!! A is not referenced.
250	!! Note that when DIAG = 'U' or 'u', the diagonal elements of
251	!! A are not referenced either, but are assumed to be unity.
252	!! Unchanged on exit.
253
254	!! LDA - INTEGER.
255	!! On entry, LDA specifies the first dimension of A as declared
256	!! in the calling (sub) program. LDA must be at least
257	!! max( 1, n ).
258	!! Unchanged on exit.
259
260	!! X - DOUBLE PRECISION array of dimension at least
261	!! ( 1 + ( n - 1 )*abs( INCX ) ).
262	!! Before entry, the incremented array X must contain the n
263	!! element right-hand side vector b. On exit, X is overwritten
264	!! with the solution vector x.
265
266	!! INCX - INTEGER.
267	!! On entry, INCX specifies the increment for the elements of
268	!! X. INCX must not be zero.
269	!! Unchanged on exit.
270
271
272	!! Level 2 Blas routine.
273
274	!! -- Written on 22-October-1986.
275	!! Jack Dongarra, Argonne National Lab.
276	!! Jeremy Du Croz, Nag Central Office.
277	!! Sven Hammarling, Nag Central Office.
278	!! Richard Hanson, Sandia National Labs.
279
280
281	!! .. Parameters ..
282	! DOUBLE PRECISION ZERO
283	! PARAMETER ( ZERO = 0.0D+0 )
284	REAL (kind=8) :: ZERO
285	PARAMETER ( ZERO = 0.0 )
286	!! .. Local Scalars ..
287	! DOUBLE PRECISION TEMP
288	REAL (kind=8) :: TEMP
289	INTEGER I, INFO, IX, J, JX, KX
290	LOGICAL NOUNIT
291	!! .. External Functions ..
292	LOGICAL LSAME
293	EXTERNAL LSAME
294	!! .. External Subroutines ..
295	EXTERNAL XERBLA
296	!! .. Intrinsic Functions ..
297	INTRINSIC MAX
298	!! ..
299	!! .. Executable Statements ..
300
301	!! Test the input parameters.
302
303	INFO = 0
304	IF ( .NOT.LSAME( UPLO , 'U' ).AND. &
305	.NOT.LSAME( UPLO , 'L' ) )THEN
306	INFO = 1
307	ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. &
308	.NOT.LSAME( TRANS, 'T' ).AND. &
309	.NOT.LSAME( TRANS, 'C' ) )THEN
310	INFO = 2
311	ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. &
312	.NOT.LSAME( DIAG , 'N' ) )THEN
313	INFO = 3
314	ELSE IF( N < 0 )THEN
315	INFO = 4
316	ELSE IF( LDA < MAX( 1, N ) )THEN
317	INFO = 6
318	ELSE IF( INCX == 0 )THEN
319	INFO = 8
320	END IF
321	IF( INFO /= 0 )THEN
322	CALL XERBLA( 'DTRSV ', INFO )
323	RETURN
324	END IF
325
326	!! Quick return if possible.
327
328	IF( N == 0 ) RETURN
329
330	NOUNIT = LSAME( DIAG, 'N' )
331
332	!! Set up the start point in X if the increment is not unity. This
333	!! will be ( N - 1 )*INCX too small for descending loops.
334
335	IF( INCX <= 0 )THEN
336	KX = 1 - ( N - 1 )*INCX
337	ELSE IF( INCX /= 1 )THEN
338	KX = 1
339	END IF
340
341	!! Start the operations. In this version the elements of A are
342	!! accessed sequentially with one pass through A.
343
344	IF( LSAME( TRANS, 'N' ) )THEN
345
346	!! Form x := inv( A )*x.
347
348	IF( LSAME( UPLO, 'U' ) )THEN
349	IF( INCX == 1 )THEN
350	DO 20, J = N, 1, -1
351	IF( X( J ) /= ZERO )THEN
352	IF( NOUNIT ) X( J ) = X( J )/A( J, J )
353	TEMP = X( J )
354	DO 10, I = J - 1, 1, -1
355	X( I ) = X( I ) - TEMP*A( I, J )
356	10 CONTINUE
357	END IF
358	20 CONTINUE
359	ELSE
360	JX = KX + ( N - 1 )*INCX
361	DO 40, J = N, 1, -1
362	IF( X( JX ) /= ZERO )THEN
363	IF( NOUNIT ) X( JX ) = X( JX )/A( J, J )
364	TEMP = X( JX )
365	IX = JX
366	DO 30, I = J - 1, 1, -1
367	IX = IX - INCX
368	X( IX ) = X( IX ) - TEMP*A( I, J )
369	30 CONTINUE
370	END IF
371	JX = JX - INCX
372	40 CONTINUE
373	END IF
374	ELSE
375	IF( INCX == 1 )THEN
376	DO 60, J = 1, N
377	IF( X( J ) /= ZERO )THEN
378	IF( NOUNIT ) X( J ) = X( J )/A( J, J )
379	TEMP = X( J )
380	DO 50, I = J + 1, N
381	X( I ) = X( I ) - TEMP*A( I, J )
382	50 CONTINUE
383	END IF
384	60 CONTINUE
385	ELSE
386	JX = KX
387	DO 80, J = 1, N
388	IF( X( JX ) /= ZERO )THEN
389	IF( NOUNIT ) X( JX ) = X( JX )/A( J, J )
390	TEMP = X( JX )
391	IX = JX
392	DO 70, I = J + 1, N
393	IX = IX + INCX
394	X( IX ) = X( IX ) - TEMP*A( I, J )
395	70 CONTINUE
396	END IF
397	JX = JX + INCX
398	80 CONTINUE
399	END IF
400	END IF
401	ELSE
402
403	!! Form x := inv( A' )*x.
404
405	IF( LSAME( UPLO, 'U' ) )THEN
406	IF( INCX == 1 )THEN
407	DO 100, J = 1, N
408	TEMP = X( J )
409	DO 90, I = 1, J - 1
410	TEMP = TEMP - A( I, J )*X( I )
411	90 CONTINUE
412	IF( NOUNIT ) TEMP = TEMP/A( J, J )
413	X( J ) = TEMP
414	100 CONTINUE
415	ELSE
416	JX = KX
417	DO 120, J = 1, N
418	TEMP = X( JX )
419	IX = KX
420	DO 110, I = 1, J - 1
421	TEMP = TEMP - A( I, J )*X( IX )
422	IX = IX + INCX
423	110 CONTINUE
424	IF( NOUNIT ) TEMP = TEMP/A( J, J )
425	X( JX ) = TEMP
426	JX = JX + INCX
427	120 CONTINUE
428	END IF
429	ELSE
430	IF( INCX == 1 )THEN
431	DO 140, J = N, 1, -1
432	TEMP = X( J )
433	DO 130, I = N, J + 1, -1
434	TEMP = TEMP - A( I, J )*X( I )
435	130 CONTINUE
436	IF( NOUNIT ) TEMP = TEMP/A( J, J )
437	X( J ) = TEMP
438	140 CONTINUE
439	ELSE
440	KX = KX + ( N - 1 )*INCX
441	JX = KX
442	DO 160, J = N, 1, -1
443	TEMP = X( JX )
444	IX = KX
445	DO 150, I = N, J + 1, -1
446	TEMP = TEMP - A( I, J )*X( IX )
447	IX = IX - INCX
448	150 CONTINUE
449	IF( NOUNIT ) TEMP = TEMP/A( J, J )
450	X( JX ) = TEMP
451	JX = JX - INCX
452	160 CONTINUE
453	END IF
454	END IF
455	END IF
456
457	END SUBROUTINE DTRSV
458	!---------------------------------------------------------
459	SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
460	!! .. Scalar Arguments ..
461	! DOUBLE PRECISION ALPHA
462	REAL (kind=8) :: ALPHA
463	INTEGER INCX, INCY, LDA, M, N
464	!! .. Array Arguments ..
465	! DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
466	REAL (kind=8) :: A( LDA, * ), X( * ), Y( * )
467	!! ..
468
469	!! Purpose
470	!! =======
471
472	!! DGER performs the rank 1 operation
473
474	!! A := alphaxy' + A,
475
476	!! where alpha is a scalar, x is an m element vector, y is an n element
477	!! vector and A is an m by n matrix.
478
479	!! Parameters
480	!! ==========
481
482	!! M - INTEGER.
483	!! On entry, M specifies the number of rows of the matrix A.
484	!! M must be at least zero.
485	!! Unchanged on exit.
486
487	!! N - INTEGER.
488	!! On entry, N specifies the number of columns of the matrix A.
489	!! N must be at least zero.
490	!! Unchanged on exit.
491
492	!! ALPHA - DOUBLE PRECISION.
493	!! On entry, ALPHA specifies the scalar alpha.
494	!! Unchanged on exit.
495
496	!! X - DOUBLE PRECISION array of dimension at least
497	!! ( 1 + ( m - 1 )*abs( INCX ) ).
498	!! Before entry, the incremented array X must contain the m
499	!! element vector x.
500	!! Unchanged on exit.
501
502	!! INCX - INTEGER.
503	!! On entry, INCX specifies the increment for the elements of
504	!! X. INCX must not be zero.
505	!! Unchanged on exit.
506
507	!! Y - DOUBLE PRECISION array of dimension at least
508	!! ( 1 + ( n - 1 )*abs( INCY ) ).
509	!! Before entry, the incremented array Y must contain the n
510	!! element vector y.
511	!! Unchanged on exit.
512
513	!! INCY - INTEGER.
514	!! On entry, INCY specifies the increment for the elements of
515	!! Y. INCY must not be zero.
516	!! Unchanged on exit.
517
518	!! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
519	!! Before entry, the leading m by n part of the array A must
520	!! contain the matrix of coefficients. On exit, A is
521	!! overwritten by the updated matrix.
522
523	!! LDA - INTEGER.
524	!! On entry, LDA specifies the first dimension of A as declared
525	!! in the calling (sub) program. LDA must be at least
526	!! max( 1, m ).
527	!! Unchanged on exit.
528
529
530	!! Level 2 Blas routine.
531
532	!! -- Written on 22-October-1986.
533	!! Jack Dongarra, Argonne National Lab.
534	!! Jeremy Du Croz, Nag Central Office.
535	!! Sven Hammarling, Nag Central Office.
536	!! Richard Hanson, Sandia National Labs.
537
538
539	!! .. Parameters ..
540	! DOUBLE PRECISION ZERO
541	! PARAMETER ( ZERO = 0.0D+0 )
542	REAL (kind=8) :: ZERO
543	PARAMETER ( ZERO = 0.0 )
544	!! .. Local Scalars ..
545	! DOUBLE PRECISION TEMP
546	REAL (kind=8) :: TEMP
547	INTEGER I, INFO, IX, J, JY, KX
548	!! .. External Subroutines ..
549	EXTERNAL XERBLA
550	!! .. Intrinsic Functions ..
551	INTRINSIC MAX
552	!! ..
553	!! .. Executable Statements ..
554
555	!! Test the input parameters.
556
557	INFO = 0
558	IF ( M < 0 )THEN
559	INFO = 1
560	ELSE IF( N < 0 )THEN
561	INFO = 2
562	ELSE IF( INCX == 0 )THEN
563	INFO = 5
564	ELSE IF( INCY == 0 )THEN
565	INFO = 7
566	ELSE IF( LDA < MAX( 1, M ) )THEN
567	INFO = 9
568	END IF
569	IF( INFO /= 0 )THEN
570	CALL XERBLA( 'DGER ', INFO )
571	RETURN
572	END IF
573
574	!! Quick return if possible.
575
576	IF( ( M == 0 ).OR.( N == 0 ).OR.( ALPHA == ZERO ) ) &
577	RETURN
578
579	!! Start the operations. In this version the elements of A are
580	!! accessed sequentially with one pass through A.
581
582	IF( INCY > 0 )THEN
583	JY = 1
584	ELSE
585	JY = 1 - ( N - 1 )*INCY
586	END IF
587	IF( INCX == 1 )THEN
588	DO 20, J = 1, N
589	IF( Y( JY ) /= ZERO )THEN
590	TEMP = ALPHA*Y( JY )
591	DO 10, I = 1, M
592	A( I, J ) = A( I, J ) + X( I )*TEMP
593	10 CONTINUE
594	END IF
595	JY = JY + INCY
596	20 CONTINUE
597	ELSE
598	IF( INCX > 0 )THEN
599	KX = 1
600	ELSE
601	KX = 1 - ( M - 1 )*INCX
602	END IF
603	DO 40, J = 1, N
604	IF( Y( JY ) /= ZERO )THEN
605	TEMP = ALPHA*Y( JY )
606	IX = KX
607	DO 30, I = 1, M
608	A( I, J ) = A( I, J ) + X( IX )*TEMP
609	IX = IX + INCX
610	30 CONTINUE
611	END IF
612	JY = JY + INCY
613	40 CONTINUE
614	END IF
615
616	END SUBROUTINE DGER
617	!---------------------------------------------------------
618	SUBROUTINE XERBLA ( SRNAME, INFO )
619	!! .. Scalar Arguments ..
620	INTEGER INFO
621	CHARACTER (len=6) :: SRNAME
622	!! ..
623
624	!! Purpose
625	!! =======
626
627	!! XERBLA is an error handler for the Level 2 BLAS routines.
628
629	!! It is called by the Level 2 BLAS routines if an input parameter is
630	!! invalid.
631
632	!! Installers should consider modifying the STOP statement in order to
633	!! call system-specific exception-handling facilities.
634
635	!! Parameters
636	!! ==========
637
638	!! SRNAME - CHARACTER*6.
639	!! On entry, SRNAME specifies the name of the routine which
640	!! called XERBLA.
641
642	!! INFO - INTEGER.
643	!! On entry, INFO specifies the position of the invalid
644	!! parameter in the parameter-list of the calling routine.
645
646
647	!! Auxiliary routine for Level 2 Blas.
648
649	!! Written on 20-July-1986.
650
651	!! .. Executable Statements ..
652
653	WRITE (*,99999) SRNAME, INFO
654
655	STOP
656
657	99999 FORMAT ( ' ** On entry to ', A6, ' parameter number ', I2, &
658	' had an illegal value' )
659
660	END SUBROUTINE XERBLA
661	!-----------------------------------------------------------
662	FUNCTION lsame( c1, c2 )
663	logical lsame
664	CHARACTER (len=*), INTENT(in) :: c1, c2
665	IF( c1 == c2 ) THEN
666	lsame=.TRUE.
667	ELSE
668	lsame=.FALSE.
669	ENDIF
670	END FUNCTION lsame

Note: See TracBrowser for help on using the repository browser.

Download in other formats: