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 @ 3

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

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

Download in other formats: