!*==sgemm.f90  processed by SPAG 7.51RB at 20:37 on  3 Mar 2022
!> \brief \b SGEMM
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!  Definition:
!  ===========
!
!       SUBROUTINE SGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
!
!       .. Scalar Arguments ..
!       REAL ALPHA,BETA
!       INTEGER K,LDA,LDB,LDC,M,N
!       CHARACTER TRANSA,TRANSB
!       ..
!       .. Array Arguments ..
!       REAL A(LDA,*),B(LDB,*),C(LDC,*)
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> SGEMM  performs one of the matrix-matrix operations
!>
!>    C := alpha*op( A )*op( B ) + beta*C,
!>
!> where  op( X ) is one of
!>
!>    op( X ) = X   or   op( X ) = X**T,
!>
!> alpha and beta are scalars, and A, B and C are matrices, with op( A )
!> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] TRANSA
!> \verbatim
!>          TRANSA is CHARACTER*1
!>           On entry, TRANSA specifies the form of op( A ) to be used in
!>           the matrix multiplication as follows:
!>
!>              TRANSA = 'N' or 'n',  op( A ) = A.
!>
!>              TRANSA = 'T' or 't',  op( A ) = A**T.
!>
!>              TRANSA = 'C' or 'c',  op( A ) = A**T.
!> \endverbatim
!>
!> \param[in] TRANSB
!> \verbatim
!>          TRANSB is CHARACTER*1
!>           On entry, TRANSB specifies the form of op( B ) to be used in
!>           the matrix multiplication as follows:
!>
!>              TRANSB = 'N' or 'n',  op( B ) = B.
!>
!>              TRANSB = 'T' or 't',  op( B ) = B**T.
!>
!>              TRANSB = 'C' or 'c',  op( B ) = B**T.
!> \endverbatim
!>
!> \param[in] M
!> \verbatim
!>          M is INTEGER
!>           On entry,  M  specifies  the number  of rows  of the  matrix
!>           op( A )  and of the  matrix  C.  M  must  be at least  zero.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>           On entry,  N  specifies the number  of columns of the matrix
!>           op( B ) and the number of columns of the matrix C. N must be
!>           at least zero.
!> \endverbatim
!>
!> \param[in] K
!> \verbatim
!>          K is INTEGER
!>           On entry,  K  specifies  the number of columns of the matrix
!>           op( A ) and the number of rows of the matrix op( B ). K must
!>           be at least  zero.
!> \endverbatim
!>
!> \param[in] ALPHA
!> \verbatim
!>          ALPHA is REAL
!>           On entry, ALPHA specifies the scalar alpha.
!> \endverbatim
!>
!> \param[in] A
!> \verbatim
!>          A is REAL array, dimension ( LDA, ka ), where ka is
!>           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
!>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
!>           part of the array  A  must contain the matrix  A,  otherwise
!>           the leading  k by m  part of the array  A  must contain  the
!>           matrix A.
!> \endverbatim
!>
!> \param[in] LDA
!> \verbatim
!>          LDA is INTEGER
!>           On entry, LDA specifies the first dimension of A as declared
!>           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
!>           LDA must be at least  max( 1, m ), otherwise  LDA must be at
!>           least  max( 1, k ).
!> \endverbatim
!>
!> \param[in] B
!> \verbatim
!>          B is REAL array, dimension ( LDB, kb ), where kb is
!>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
!>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
!>           part of the array  B  must contain the matrix  B,  otherwise
!>           the leading  n by k  part of the array  B  must contain  the
!>           matrix B.
!> \endverbatim
!>
!> \param[in] LDB
!> \verbatim
!>          LDB is INTEGER
!>           On entry, LDB specifies the first dimension of B as declared
!>           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
!>           LDB must be at least  max( 1, k ), otherwise  LDB must be at
!>           least  max( 1, n ).
!> \endverbatim
!>
!> \param[in] BETA
!> \verbatim
!>          BETA is REAL
!>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
!>           supplied as zero then C need not be set on input.
!> \endverbatim
!>
!> \param[in,out] C
!> \verbatim
!>          C is REAL array, dimension ( LDC, N )
!>           Before entry, the leading  m by n  part of the array  C must
!>           contain the matrix  C,  except when  beta  is zero, in which
!>           case C need not be set on entry.
!>           On exit, the array  C  is overwritten by the  m by n  matrix
!>           ( alpha*op( A )*op( B ) + beta*C ).
!> \endverbatim
!>
!> \param[in] LDC
!> \verbatim
!>          LDC is INTEGER
!>           On entry, LDC specifies the first dimension of C as declared
!>           in  the  calling  (sub)  program.   LDC  must  be  at  least
!>           max( 1, m ).
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup single_blas_level3
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>  Level 3 Blas routine.
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE SGEMM(Transa,Transb,M,N,K,Alpha,A,Lda,B,Ldb,Beta,C,Ldc)
      IMPLICIT NONE
!*--SGEMM191
!
!  -- Reference BLAS level3 routine (version 3.7.0) --
!  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      REAL Alpha , Beta
      INTEGER K , Lda , Ldb , Ldc , M , N
      CHARACTER Transa , Transb
!     ..
!     .. Array Arguments ..
      REAL A(Lda,*) , B(Ldb,*) , C(Ldc,*)
!     ..
!
!  =====================================================================
!
!     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL LSAME
!     ..
!     .. External Subroutines ..
      EXTERNAL XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MAX
!     ..
!     .. Local Scalars ..
      REAL temp
      INTEGER i , info , j , l , nrowa , nrowb
      LOGICAL nota , notb
!     ..
!     .. Parameters ..
      REAL ONE , ZERO
      PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
!     ..
!
!     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
!     transposed and set  NROWA and NROWB  as the number of rows of  A
!     and  B  respectively.
!
      nota = LSAME(Transa,'N')
      notb = LSAME(Transb,'N')
      IF ( nota ) THEN
         nrowa = M
      ELSE
         nrowa = K
      ENDIF
      IF ( notb ) THEN
         nrowb = K
      ELSE
         nrowb = N
      ENDIF
!
!     Test the input parameters.
!
      info = 0
      IF ( (.NOT.nota) .AND. (.NOT.LSAME(Transa,'C')) .AND.             &
     &     (.NOT.LSAME(Transa,'T')) ) THEN
         info = 1
      ELSEIF ( (.NOT.notb) .AND. (.NOT.LSAME(Transb,'C')) .AND.         &
     &         (.NOT.LSAME(Transb,'T')) ) THEN
         info = 2
      ELSEIF ( M<0 ) THEN
         info = 3
      ELSEIF ( N<0 ) THEN
         info = 4
      ELSEIF ( K<0 ) THEN
         info = 5
      ELSEIF ( Lda<MAX(1,nrowa) ) THEN
         info = 8
      ELSEIF ( Ldb<MAX(1,nrowb) ) THEN
         info = 10
      ELSEIF ( Ldc<MAX(1,M) ) THEN
         info = 13
      ENDIF
      IF ( info/=0 ) THEN
         CALL XERBLA('SGEMM ',info)
         RETURN
      ENDIF
!
!     Quick return if possible.
!
      IF ( (M==0) .OR. (N==0) .OR.                                      &
     &     (((Alpha==ZERO) .OR. (K==0)) .AND. (Beta==ONE)) ) RETURN
!
!     And if  alpha.eq.zero.
!
      IF ( Alpha==ZERO ) THEN
         IF ( Beta==ZERO ) THEN
            DO j = 1 , N
               DO i = 1 , M
                  C(i,j) = ZERO
               ENDDO
            ENDDO
         ELSE
            DO j = 1 , N
               DO i = 1 , M
                  C(i,j) = Beta*C(i,j)
               ENDDO
            ENDDO
         ENDIF
         RETURN
      ENDIF
!
!     Start the operations.
!
      IF ( notb ) THEN
         IF ( nota ) THEN
!
!           Form  C := alpha*A*B + beta*C.
!
            DO j = 1 , N
               IF ( Beta==ZERO ) THEN
                  DO i = 1 , M
                     C(i,j) = ZERO
                  ENDDO
               ELSEIF ( Beta/=ONE ) THEN
                  DO i = 1 , M
                     C(i,j) = Beta*C(i,j)
                  ENDDO
               ENDIF
               DO l = 1 , K
                  temp = Alpha*B(l,j)
                  DO i = 1 , M
                     C(i,j) = C(i,j) + temp*A(i,l)
                  ENDDO
               ENDDO
            ENDDO
         ELSE
!
!           Form  C := alpha*A**T*B + beta*C
!
            DO j = 1 , N
               DO i = 1 , M
                  temp = ZERO
                  DO l = 1 , K
                     temp = temp + A(l,i)*B(l,j)
                  ENDDO
                  IF ( Beta==ZERO ) THEN
                     C(i,j) = Alpha*temp
                  ELSE
                     C(i,j) = Alpha*temp + Beta*C(i,j)
                  ENDIF
               ENDDO
            ENDDO
         ENDIF
      ELSEIF ( nota ) THEN
!
!           Form  C := alpha*A*B**T + beta*C
!
         DO j = 1 , N
            IF ( Beta==ZERO ) THEN
               DO i = 1 , M
                  C(i,j) = ZERO
               ENDDO
            ELSEIF ( Beta/=ONE ) THEN
               DO i = 1 , M
                  C(i,j) = Beta*C(i,j)
               ENDDO
            ENDIF
            DO l = 1 , K
               temp = Alpha*B(j,l)
               DO i = 1 , M
                  C(i,j) = C(i,j) + temp*A(i,l)
               ENDDO
            ENDDO
         ENDDO
      ELSE
!
!           Form  C := alpha*A**T*B**T + beta*C
!
         DO j = 1 , N
            DO i = 1 , M
               temp = ZERO
               DO l = 1 , K
                  temp = temp + A(l,i)*B(j,l)
               ENDDO
               IF ( Beta==ZERO ) THEN
                  C(i,j) = Alpha*temp
               ELSE
                  C(i,j) = Alpha*temp + Beta*C(i,j)
               ENDIF
            ENDDO
         ENDDO
      ENDIF
!
!
!     End of SGEMM .
!
      END SUBROUTINE SGEMM


text:
Shahs Halal Food is a Restaurants food trailer on the UC Davis campus that serves sandwiches and rice platters with grilled lamb or chicken, similar in style to the halal food carts of NYC such as the wiki:Wikipedia:53rd and 6th famous cart on 53rd and 6th. Prices for rice platters include a soft drink. There are three choices of sauce, http://dinersjournal.blogs.nytimes.com/2009/06/02/whatsinthewhitesaucehedoesntwanttoknow/ white sauce, hot sauce (Sriracha), and BBQ sauce.

20091001 13:51:49 nbsp Tried the mixed (lamb and chicken) rice platter today.  Overall, quite good and similar in quality to the ones in New York City.  You get white sauce automatically and will add Sriracha as hot sauce (they dont seem to make their own).  I hope they make it and hope to see more of these types of campus food options. Users/spencerchang

20091005 16:06:38 nbsp is this actually a campus sanctioned food vendor or just some random guy selling food on campus? Users/StevenLee

20091006 22:03:01 nbsp They ended up forgetting my order and two people in front of me as well. I waited around 10 extra minutes for the lamb gyro sandwich. The wait was worth it though, everything was hot and very flavorful. Not too expensive either. Users/KevinAdamski

20091007 00:46:58 nbsp How does this place compare to Sams and Ali Baba?

I saw that they include a drink in the price.  Is it possible to get food from them for less without the drink?  (Not everyone wants to drink soda or have a cup of liquid to carry around.) Users/IDoNotExist

20091008 21:25:17 nbsp I saw this place a few weeks back and I was really pumped to try it. Alas, it sucked. I should have just walked away when the dude taking orders apparently thinks I want a jyerow. I dont mean to be a jerk, but Im sorry, if what you do is sell gyros please AT LEAST recognize the correct pronunciation is NOT jyerows.  The meat wasnt bad but the gyros was just meat, lettuce and tzatziki.  For real is that it?  Not the most intimidating looking gyros Ive ever seen either. Oh well. Users/Patrick.M

20091009 13:09:01 nbsp There are several different correct pronunciations for gyros, including yeeros, jierow, and geeros (with a hard g). Users/IDoNotExist

20091015 16:54:50 nbsp Bring a bathing suit if you eat at this place, because the lamb, while tasty, is literally swimming in yogurt sauce.  The food wasnt bad, although it doesnt compare to Sams in quality or taste.  The gyro is more like what you would get at Sams than at Ali Baba, but Sams is definitely better.  Shahs lamb is rather greasy.  While I liked the taste of what I got there, I could feel the grease spreading down through my stomach, and I felt kinda greasy on the inside hours after eating there.  If you were going to have to choose between Shahs and the Silo, which is a short walk away, Id have a hard time recommending one over the other.  Both are a bit overpriced and not very healthy.  But if you werent restricted to eating in that general area, Id definitely suggest taking a few minutes to bike to Sams or Ali Baba instead.

On the IDoNotExist SushiMeter (TM), this place scores rather poorly, as it doesnt serve sushi.  On the IDoNotExist MediterraneanFoodMeter (TM), this place scores a 4/10 (tasty food, but unhealthy, WAY too much yogurt sauce, and unpleasant greasy aftereffects...) Users/IDoNotExist

20091017 10:24:19 nbsp I Loved their food it was so good. Much better than Sams. Users/CRITIC

20091019 09:42:41 nbsp I went to this spot on the corner of Hutchinson and California for the first time on Friday. I got the mixed rice plate and it was great. Also, a lot of food for the $6.50 price tag. I stood in line for about 3 mins and got my meal in 2 more, service was fast. UMMMMMM, was the white sauce good; a little heavy but very good.  Monday, ima try the Gyro, Eruo or however it is pronounced.  Users/Pinnacle

20091206 21:04:19 nbsp cant really compare this place to a restaurant, its a foodcart afterall. id say the lamb gyro was pretty good considering it was premade and served to me in less than 5 mins. price wasnt that bad either considering all the overpriced dishes here in davis Users/betabeta

20091226 13:56:44 nbsp Someones deity clearly blessed me for finally providing something that isnt the Silo to eat. Im in the Physics building and dont entertain the idea of going the MU or off campus to eat usually. Shahs is new, extremely tasty, and I do not foresee getting sick of it any time soon. Ive eaten it 3 times in a work week and been excited every time. You can eat for cheaper at Taco Bell if you wanted... but why? If you have tried this place, do so. The lamb is processed, but good. The chicken  happy me. Users/MichaelAWoods

20091227 14:27:26 nbsp I for one love the yogurt sauce, the guys who run this trailer are pretty cool! The food is amazing, I love getting all my orders spicy. They chicken gyro is amazing, and the lamb rice platter was pretty solid as well. Users/StevenDaubert

20100113 15:52:33 nbsp I came here today and got the lamb over rice.  Was a bit pricey but I was satisfied.  Maybe I am not familiar enough with this kind of food, but the lamb didnt taste or feel very lamby to me.  Despite that the meal was very tasty and I will most likely return to try the chicken over rice.

I asked them not to put in the hot sauce, but they accidentally put it anyways and were about to remake the meal before I told them its okay.  I was glad the sauce was there, because I think the meal wouldve been dry otherwise.  If you dont mind the spice and calories, make sure to get the sauces. Users/BrandonWong

20100203 10:45:45 nbsp I liked it.  Lamb over rice was great, and is enough to be filling when I want a meal.  The orders (or at least mine) come with a drink, too.  All food is wrapped to go, of course.  It was reasonably fast, too.  Will definitely be going there again.   Users/KayChang

20100529 15:01:08 nbsp I really like this place although Ive only tried the lamb over rice and chicken over rice. I think the portion is really generous. For lunch, my boyfriend and I can share one and still be full. And for dinner, I can only eat like 2/3 of it while my boyfriend can barely finish it. The bad part is the onions in the salad which I cant ask for them to leave it out. Users/snsdgee

20100529 17:48:26 nbsp Chicken Gyro is still very bomb, highly recommended Users/StevenDaubert

20100608 16:55:03 nbsp The chicken gyro is the best thing you can order on campus. so tasty!! Users/KatieP

20100630 16:49:47 nbsp Messy, spicy gyro goodness. Great deal and they are always super nice.  Users/LokiAbbi

20100706 20:29:33 nbsp Can someone please tell me where they moved?? The school construction has now expanded to that corner. Users/WengFF

20100707 15:14:32 nbsp They moved to the grassy area in front of Hutchison Hall on California Ave.  Its just south of Admundson Hall and across the street from Kerr Hall. Users/JoeCastillo

20100714 13:02:22 nbsp I keep hearing rumors that they have falafel.  Has anyone tried it if so? Users/IAmNewHere

20100715 12:31:39 nbsp Answered my own question about falafel: they have it and its good.  I still prefer Sams, which gives more for the same price (though drink is included here and extra at Sams), and I prefer Sams flavor just a bit more.  Definitely would come here over Ali Baba or Café Méditerranée for falafel, cant really speak to anything else.   Users/IAmNewHere

20100802 15:56:15 nbsp Overall, this cart is great. Quite delicious. Today I got a chicken YEERO ( ;D ) although what I really wanted was a lamb one. I asked the guy if it was really lamb or lambbeef mixture, alas it was the latter (read labels, folks!). This isnt their fault, Kronos sends out the pillars of gyro meat from the factory that way, but I hate beef so opted for the slightly less reprehensible (food, inc. anyone?) chicken. Also, while the food may be halaal it is definitely not kosher, if that is a concern for you. But to return to my first point, Ill definitely be hitting this place up again, at least until 2019 when they finish revamping the coffee mall.  Users/YesItsMe

20100902 11:54:39 nbsp does anyone know where theyre located now? Users/MBui
  South East side of the Silo Daubert

20100929 22:50:44 nbsp I had the chicken gyro and chicken over rice here last year. They were delicious. Yes, the gyro is a bit messy and greasy but SO good!  Today I had the chicken over rice again, and it seemed to me like the yogurt (or white sauce?) was a bit tangier than last year? At any rate, I always ask for hot sauce and it comes standard with the yogurt (or again, is it yogurt or white sauce?). I definitely recommend! Users/redbike

20100930 09:29:21 nbsp I could eat this food every day! I go to the stand on a pretty regular basis (so much so that the cashier guy recognizes me and kind of laughs when I tell him what I wantthe same thing every time). Ive tried everything on the menu, but my favourite is definitely the Lamb over rice! The chicken is not excellent in my opinion, its just like any chicken youd find in a cheap taco joint. However the lamb is delicious! Most people complain that the lamb is ground up and not whole pieces, but I prefer it this wasI mean come on, would you really want whole chunks of lamb in there? Theyd be tough, less flavorful, and cumbersome to eat.

Anyway, Id recommend this place a million times over, and I think that the foodis definitely worth the price (hey, its usually at least 2 servings worth for me)! Users/AshleyMatson

20100930 22:27:23 nbsp I really love this place. Ive only ever ordered the chicken over rice, and I think its delicious, especially with extra spicy sauce. Ive wanted to try the gyros, but I always end up falling for the platter. And for under seven dollars, the dish serves up a main course (the chicken), a side (the rice), a salad of sorts, and a  soda, so not too bad, but a little more than I can afford on a daily basis. Plus the guys that work there are really nice, and you get your food quickly even when the line is long. Users/jwieland1989

20101020 00:39:25 nbsp The majority of people get the chicken over rice (I know this because I was waiting in their long line usually gets long on the hour and listening to what other people ordered) Out of ten people about 7 got the chicken over rice. IMO I think its the best dish on the menu.  If youre a vegetarian, they can accommodate you as well.
I am addicted to Sri Racha Red Rooster Sauce and I love the fact that they (upon request...ie Can I get extra hot sauce?) slather your dish with it.
The white sauce is not spectacular (and rather greasy), I once peaked inside the truck and saw hidden valley ranch bottles so Im guessing they either 1)only use hidden valley ranch as their sole white sauce or
2)they mix the ranch with yogurt.
The meat is tender and flavorful (except for Lamb).  The rice is sometimes a little too firm (I prefer my rice to be more fluffy).
Overall I go here because of the convenience and because of the customer service.  The price is a little steep but for the quality of the dish its well worth it.

Now, I wish we had more truck foods on campus...I believe there is a petition on FB to get a taco truck.
The link is :http://www.facebook.com/home.php?skgroup_162048093818919
If it doesnt work, search Bring a Taco Truck to UC Davis!! in the FB search bar.
As much as I love Shahs, I need to have a diff. and cheaper place to go in between visits to Shahs. Users/rootmeansquare
  This would be so much easier if we could loosen Sodexos grip on us. I seriously do not understand why people have to pay them to operate on campus. Utter nonsense. The higherups who allow this need to be force fed the garbage they serve in the dining commons for over ten dollars a meal. Users/hankim

20111114 20:35:14 nbsp I usually get the Lamb over fanta NO SPICYYYYY!!! YUM!!! This place is good but sometimes Travis shows up Users/KCSlater

20120301 13:59:21 nbsp In my opinion, this is the best food on campus by a long shot. I eat here nearly everyday and always get plenty of great food with excellent service. The cost is comparable to all other campus food venues considering how much you get, and a drink is always included. I recommend their food to everyone! Keep up the great work, Shahs!  Users/JohnSchreiber

20120301 17:03:43 nbsp The chicken gyro, spicy, always makes my day. Om NOM.

I always look down upon those eating from the fast food places at the Silo whenever I get food from Shahs. I just wish there was a taco truck on campus as well. I know it would do killer business.  Users/CecilioPadilla

20120620 13:58:20 nbsp Very good food, I agree with most of the comments.  Sadly, I will not be going back.  Today I bit into a piece of glass I found in my food, then cut my finger removing the glass from my mouth.  I took it back and got a refund, but they pretty much just shrugged it off, saying, it happens..  They should really be more careful with what theyre feeding people.  Im really disappointed.  This was my favorite campus food. Users/NancyPants

20131203 23:38:10 nbsp I love this place, its sooooo good. Prices are reasonable but buying too often will be pricey. I always get the chicken over rice with extra meat. Gets me SO FULL. The wait is quick they literally give you your food in about a minute of ordering and it isnt cold. If I knew I wouldnt get sick of this place I would eat there everyday! This is a must try at least once if you are in Davis. Users/LILP

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE RD_SYNT ( gemftm, jf6, jftmst, jftmen, rdata,
     +			     gemstm, ps12tb, iret )
C************************************************************************
C* RD_SYNT								*
C*									*
C* This subroutine converts the GEMPAK date/time array into synoptic	*
C* times.  It also sets up the position moving table for 12-hour	*
C* parameters.								*
C*									*
C* RD_SYNT  ( GEMFTM, JF6, JFTMST, JFTMEN, RDATA, GEMSTM, PS12TB, IRET )*
C*									*
C* Input parameters:							*
C*      GEMFTM (*)	CHAR*		GEMPAK date/time		*
C*      JF6		INTEGER		Index of 1st hour in extended pt*
C*      JFTMST		INTEGER		Index of 1st valid GEMFTM	*
C*									*
C* Input and output parameters:						*
C*      JFTMEN		INTEGER		Index of last GEMFTM		*
C*      RDATA		REAL		Forecast data			*
C*	   ( MMPARM, * )						*
C*									*
C* Output parameters:							*
C*      GEMSTM (*)	CHAR*		GEMPAK Synoptic date/time	*
C*      PS12TB (*)	INTEGER		Position moving table for 12-hr *
C*	IRET		INTEGER		Return code			*
C*					   0 = normal return		*
C*					  -1 = Bad fcst hrs in 2nd part *
C**									*
C* Log:									*
C* F. J. Yen/NCEP	 9/02						*
C* F. J. Yen/NCEP	11/02	Use actual time, not 3-hr syn for table;*
C*                          	handled overflow for 12-hr parameter.	*
C************************************************************************
	INCLUDE         'GEMPRM.PRM'
C*
	CHARACTER*(*)	gemftm(*), gemstm (*)
	INTEGER		ps12tb (*)
	REAL		rdata ( MMPARM, * )
C*
	CHARACTER	gemttm*15, gem6tm*15
	INTEGER		idtarr (5), idtars (5), idfovf (5)
C------------------------------------------------------------------------
	iret = 0
C
C*	Initialize position moving table for 12-hour parameters
C
	DO ii = 1, jftmen + 1
	    ps12tb ( ii ) = 0
	END DO
C
C*	Convert the GEMPAK date/time array into synoptic times.
C*	First convert the first part to 3-hour synoptic times.
C*	Also set up the 12 hour position moving table ps12tb.
C
	IF ( jf6 .eq. 0 ) THEN
	    j1 = jftmen
	  ELSE
	    j1 = jf6 - 1
	END IF
	DO ii = jftmst, j1
C
C*	    Convert to 3-hour synoptic times
C
	    CALL TI_CTOI ( gemftm (ii), idtarr, ier )
	    DO j = 1, 5
		idtars ( j ) = idtarr ( j )
	    END DO
	    irm3 = MOD ( idtarr (4), 3)
	    IF ( irm3 .eq. 1 ) THEN
		CALL TI_SUBM ( idtarr, 60, idtarr, ier )
	      ELSE IF ( irm3 .eq. 2 ) THEN
		CALL TI_ADDM ( idtarr, 60, idtarr, ier )
	    END IF
	    CALL TI_ITOC ( idtarr, gemstm ( ii ), ier)
C
C*	    Set up position moving table
C
	    irm6 = MOD ( idtars (4), 6)
	    IF ( irm6 .eq. 2) THEN
		ps12tb ( ii ) = -1
	      ELSE IF ( irm6 .eq. 3 .or. irm6 .eq. 4 ) THEN
		ps12tb ( ii ) = 1
	    END IF
	END DO
	jf6n = jf6
	jf6c = jf6
	IF ( jf6 .ne. 0 .and. jf6. le. jftmen ) THEN
C
C*	    Have extended part; convert to 6-hour synoptic times.
C*	    Convert the first fcst hour in the extended section.
C
	    CALL TI_CTOI ( gemftm (jf6), idtarr, ier )
	    irm6 = MOD ( idtarr (4), 6 )
	    IF ( irm6 .ge. 3 ) THEN
		im = (6 - irm6) * 60
		CALL TI_ADDM ( idtarr, im, idtarr, ier)
	      ELSE IF ( irm6 .ne. 0 ) THEN
		im = irm6 * 60
		CALL TI_SUBM ( idtarr, im, idtarr, ier)
	    END IF
	    CALL TI_ITOC ( idtarr, gemttm, ier )
C
C*	    Check if have overflow.  There is overflow when
C*	    the 12-hour parameters in the last forecast hour
C*	    of the first section are to be moved to the next
C*	    (non-existent) 3-hour forecast hour.
C
	    IF ( ps12tb ( j1 ) .eq. 1 ) THEN
C
C*		There is overflow; determine the 6 hour
C*		synoptic time of the overflow and
C*		compare with the first 6 hour synoptic
C*		time of the extended section.
C
		CALL TI_CTOI ( gemftm (j1), idfovf, ier )
	        irm6 = MOD ( idfovf (4), 6 )
	        IF ( irm6 .ge. 3 ) THEN
		    im = (6 - irm6) * 60
		    CALL TI_ADDM ( idfovf, im, idfovf, ier)
	          ELSE IF ( irm6 .ne. 0 ) THEN
		    im = irm6 * 60
		    CALL TI_SUBM ( idfovf, im, idfovf, ier)
	        END IF
		CALL TI_ITOC ( idfovf, gem6tm, ier )
		CALL TI_DIFF ( gemttm, gem6tm, mindif, ier)
		IF ( mindif .gt. 0 ) THEN
C
C*		    Make room for the overflow
C
		    gemstm ( jf6n ) = gem6tm
		    jf6n = jf6n + 1
		    gemstm ( jf6n ) = gemttm
		    jftmen = jftmen + 1
		    DO mm = jftmen, jf6n, -1
			DO mp = 1, MMPARM
			    rdata ( mp, mm ) = rdata ( mp, mm-1 )
			END DO
		    END DO
		    DO mp = 1, MMPARM
			rdata ( mp, mm ) = RMISSD
		    END DO
		  ELSE IF ( mindif .eq. 0 ) THEN
		    gemstm (jf6n) = gemttm
		  ELSE
C
C*		    Error in forecast time for second section
C
		    iret = -1
		    RETURN
		END IF
	      ELSE
		    gemstm (jf6n) = gemttm
	    END IF

	    IF ( jf6c .lt. jftmen ) THEN
	        jf6n = jf6n + 1
	        DO ii = jf6n, jftmen
		    jf6c = jf6c + 1
		    CALL TI_CTOI ( gemftm (jf6c), idtarr, ier )
		    irm6 = MOD ( idtarr (4), 6 )
		    IF ( irm6 .ge. 3 ) THEN
		        im = (6 - irm6) * 60
		        CALL TI_ADDM ( idtarr, im, idtarr, ier)
		      ELSE IF ( irm6 .ne. 0 ) THEN
		        im = irm6 * 60
		        CALL TI_SUBM ( idtarr, im, idtarr, ier)
		    END IF
		    CALL TI_ITOC ( idtarr, gemstm (ii), ier)
	        END DO
	    END IF
	  ELSE IF ( ps12tb ( j1 ) .eq. 1 ) THEN
C
C*	    Overflow, but no extended section
C
            gemstm ( jf6 ) = gem6tm
            jftmen = jftmen + 1
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
#version 120

varying vec2 p;

uniform float param = 128.;

int mandel(vec2 c, int n)
{
    float res_r = 0.0;
    float res_i = 0.0;
    float res_r_tmp;
    int i = 0;
    while ((res_r * res_r + res_i * res_i < 4) && (i < n))
    {
        res_r_tmp = res_r;
        res_r = res_r * res_r - res_i * res_i + c.x;
        res_i = 2.0 * res_r_tmp * res_i + c.y;
        ++i;
    }
    return i;
}

#define KAHAN_ITER(I) \
    { \
        float y = I - c_; \
        float t = sum_ + y; \
        c_ = (t - sum_) - y; \
        sum_ = t; \
    }

float kahan_sum6(float i0, float i1, float i2, float i3, float i4, float i5)
{
    float sum_ = 0.0;
    float c_ = 0.0;
    KAHAN_ITER(i0)
    KAHAN_ITER(i1)
    KAHAN_ITER(i2)
    KAHAN_ITER(i3)
    KAHAN_ITER(i4)
    KAHAN_ITER(i5)
    return sum_;
}

#define KAHAN_SUM7(O1, O2, I0, I1, I2, I3, I4, I5, I6)\
{\
    float sum_ = 0.0;\
    float c_ = 0.0;\
    KAHAN_ITER(I0)\
    KAHAN_ITER(I1)\
    KAHAN_ITER(I2)\
    KAHAN_ITER(I3)\
    KAHAN_ITER(I4)\
    KAHAN_ITER(I5)\
    KAHAN_ITER(I6)\
    O1 = sum_; O2 = c_;\
}

#define KAHAN_SUM5(O1, O2, I0, I1, I2, I3, I4)\
{\
    float sum_ = 0.0;\
    float c_ = 0.0;\
    KAHAN_ITER(I0)\
    KAHAN_ITER(I1)\
    KAHAN_ITER(I2)\
    KAHAN_ITER(I3)\
    KAHAN_ITER(I4)\
    O1 = sum_; O2 = c_;\
}

#define SUM_SQR(A,B) A*A, 2*A*B, B*B

int mandel2(vec2 c, int n)
{
    float res_r0 = 0.0, res_r1 = 0.0;
    float res_i0 = 0.0, res_i1 = 0.0;

    float res_r0_tmp, res_r1_tmp;
    int i = 0;
    while ((kahan_sum6(SUM_SQR(res_r0,res_r1), SUM_SQR(res_i0,res_i1)) < 4) && (i < n))
    {
        res_r0_tmp = res_r0;
        res_r1_tmp = res_r1;

        KAHAN_SUM7(res_r0, res_r1,
            res_r0*res_r0, 2*res_r0*res_r1, res_r1*res_r1,
            -res_i0*res_i0, -2*res_i0*res_i1, -res_i1*res_i1,
            c.x);

        KAHAN_SUM5(res_i0, res_i1,
            2.0*res_r0_tmp*res_i0, 2.0*res_r1_tmp*res_i0, 2.0*res_r0_tmp*res_i1, 2.0*res_r1_tmp*res_i1,
            c.y);

        ++i;
    }
    return i;
}

void main()
{
    vec2 pp = p;
    pp.x = p.x * 1.5 - 0.5;
    pp.y *= 1.5;
    float color1 = 1. - mandel(pp, int(param)) / param;
    float color2 = 1. - mandel2(pp, int(param)) / param;
    gl_FragColor = vec4(color1, color2, 0.0, 1.);
}

max_stars_repo_name: hedlx/machuchu

max_issues_repo_name: hedlx/machuchu

max_forks_repo_name: hedlx/machuchu


text:
      SUBROUTINE PHO_SASDIR(X,Q2,P2,Q02,Xpga)
      IMPLICIT NONE
      REAL aem2pi , cgam , P2 , pmb , pmc , Q02 , Q2 , X , Xpga , xtmp
      INTEGER kf , kfl
C...Purpose: to evaluate the direct contribution, i.e. the C^gamma term,
C...as needed in MSbar parametrizations.
      SAVE
      DIMENSION Xpga(-6:6)
      DATA pmc/1.3/ , pmb/4.6/ , aem2pi/0.0011614/

C...Reset output.
      DO kfl = -6 , 6
         Xpga(kfl) = 0.
      END DO

C...Evaluate common x-dependent expression.
      xtmp = (X**2+(1.-X)**2)*(-LOG(X)) - 1.
      cgam = 3.*aem2pi*X*(xtmp*(1.+P2/(P2+Q02))+6.*X*(1.-X))

C...d, u, s part by simple charge factor.
      Xpga(1) = (1./9.)*cgam
      Xpga(2) = (4./9.)*cgam
      Xpga(3) = (1./9.)*cgam

C...Also fill for antiquarks.
      DO kf = 1 , 5
         Xpga(-kf) = Xpga(kf)
      END DO

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


! { dg-do run }
! { dg-options "-O0" }
!
! Test fix for PR18022.
!
! Contributed by Paul Thomas <pault@gcc.gnu.org>
!
program assign_func_dtcomp
 implicit none
 type                         ::  mytype
   real                       ::  x
   real                       ::  y
 end type mytype
 type (mytype), dimension (4) ::  z

 type                         ::  thytype
   real                       ::  x(4)
 end type thytype
 type (thytype)               ::  w
 real, dimension (4)          ::  a = (/1.,2.,3.,4./)
 real, dimension (4)          ::  b = (/5.,6.,7.,8./)

! Test the original problem is fixed.
 z(:)%x = foo (a)
 z(:)%y = foo (b)

 if (any(z%x.ne.a).or.any(z%y.ne.b)) STOP 1

! Make sure we did not break anything on the way.
 w%x(:) = foo (b)
 a = foo (b)

 if (any(w%x.ne.b).or.any(a.ne.b)) STOP 2

contains

 function foo (v) result (ans)
   real, dimension (:), intent(in)   ::  v
   real, dimension (size(v))  ::  ans
   ans = v
 end function foo

end program assign_func_dtcomp


#if defined(__CUDA)
module cholesky
use cudafor
use cublas
use cusolverDn

logical, save         :: cusolver_init=.false.
type(cusolverDnHandle):: handle_cusolver
integer, device       :: devInfo_d
integer               :: istat
complex(8), parameter::  cone=cmplx( 1.0d+0, 0.0d+0 ), czero=cmplx( 0.0D+0, 0.0D+0 )
real(8)   , parameter::  one=1.0d+0, zero=0.0D+0

! Generic interfaces for compute:
!  Cholesky factorization:
!     dpotrf(double precision)  and zpotrf( double complex matrix)
!  Inverse of a triangular matrix:
!    dtrtri(double precision)  and ztrtri( double complex matrix)
! For the TRTRI routines, only the path used by Quantum-Espresso ( Lower, non-unit)
! is implemented on GPU.
! The CPU routines will be in the BLAS library, we are just specifying the interface
! to dispatch to the CPU or GPU based on the properties of the matrix A

interface dpotrf
    subroutine dpotrf( uplo, n, A, lda, info )
      use kinds, only : dp
      character:: uplo
      integer:: n, lda,info
      real(dp):: A(lda,*)
    end subroutine dpotrf

    module procedure dpotrf_gpu
end interface dpotrf

interface zpotrf
    subroutine zpotrf( uplo, n, A, lda, info )
      use kinds, only : dp
      character:: uplo
      integer :: n, lda,info
      complex(dp):: A(lda,*)
    end subroutine zpotrf

    module procedure zpotrf_gpu
end interface zpotrf

interface dtrtri
    subroutine dtrtri( uplo, diag, n, A, lda, info )
      use kinds, only : dp
      character:: uplo, diag
      integer:: n, lda, info
      real(dp):: A(lda,*)
    end subroutine dtrtri

    module procedure dtrtri_gpu
end interface dtrtri

interface ztrtri
    subroutine ztrtri( uplo, diag, n, A, lda, info )
      use kinds, only : dp
      character:: uplo, diag
      integer :: n, lda,info
      complex(dp):: A(lda,*)
    end subroutine ztrtri

    module procedure ztrtri_gpu
end interface ztrtri

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE dpotrf_gpu( uplo, n, A, lda, info )
      use kinds, only : dp
      implicit none
      character::          UPLO
      integer :: n,lda,info
      real(dp), intent(inout), device:: A(lda,*)
      !
      ! given a matrix A, returns the Cholesky decomposition of A
      ! for cholesky matrices
      !
      real(dp),  device:: work(1)
      integer :: lwork=1
      logical:: upper
      logical,external:: lsame

      if (.not. cusolver_init) then
        istat = cusolverDnCreate(handle_cusolver)
        if (istat /= CUSOLVER_STATUS_SUCCESS) write(*,*) 'cusolver handle creation failed'
      end if

      info = -1
      upper = lsame( uplo, 'U' )
      if(upper) then
       istat= cusolverDnDpotrf(handle_cusolver, CUBLAS_FILL_MODE_UPPER, n, A, lda, work, lwork, devInfo_d)
      else
       istat= cusolverDnDpotrf(handle_cusolver, CUBLAS_FILL_MODE_LOWER, n, A, lda, work, lwork, devInfo_d)
      endif
      info=devInfo_d

    END SUBROUTINE dpotrf_gpu

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     SUBROUTINE zpotrf_gpu( uplo, n, A, lda, info )
      use kinds, only : dp
      implicit none
      character::          UPLO
      integer :: n,lda,info
      complex(dp), intent(inout), device:: A(lda,*)
      !
      ! given a matrix A, returns the Cholesky decomposition of A
      ! for cholesky matrices
      !
      complex(dp),  device:: work(1)
      integer :: lwork=1
      logical:: upper
      logical,external:: lsame

      if (.not. cusolver_init) then
        istat = cusolverDnCreate(handle_cusolver)
        if (istat /= CUSOLVER_STATUS_SUCCESS) write(*,*) 'cusolver handle creation failed'
      end if

      info = -1
      upper = lsame( uplo, 'U' )
      if(upper) then
       istat= cusolverDnZpotrf(handle_cusolver, CUBLAS_FILL_MODE_UPPER, n, A, lda, work, lwork, devInfo_d)
      else
       istat= cusolverDnZpotrf(handle_cusolver, CUBLAS_FILL_MODE_LOWER, n, A, lda, work, lwork, devInfo_d)
      endif
      info=devInfo_d

    END SUBROUTINE zpotrf_gpu

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! The two following routines are the GPU implementation of the ZTRTRI LAPACK routine.
! This routine computes the inverse of a complex lower non-unit triangular matrix.
! This is the unblocked code that operates on a block 32x32
   attributes(global) subroutine ztrti2_gpu(n,a,lda)
      use cudafor
      implicit none
      integer, value    :: n,lda
      complex(8),device :: a(lda,*)
      complex(8),shared :: a_s(32,32)
      complex(8)        :: cv
      integer           :: tx,ty,tl,i,j,ii

     if (n == 0 ) return

      tx=threadIdx%x
      ty=threadIdx%y
      ! Linear id of the thread (tx,ty)
      tl=tx+ blockDim%x*(ty-1)

      ! Load a_d in shared memory
      if (tx <= n .and. ty <= n) then
         a_s(tx,ty)=a(tx,ty)
      endif

      call syncthreads()

     ! Compute all the diagonal elements

     if (tl <=n ) a_s(tl,tl)=cone/a_s(tl,tl)

      call syncthreads()

     ! For each column working backward
      do i=n-1,1,-1

         if ( tl >i .and. tl <= n) then
            cv=czero
             do j= i+1,tl
               cv = cv + a_s(j , i ) * a_s(tl ,j)
             end do
         end if
!        call syncthreads()
         if ( tl >i .and. tl <= n) then
           a_s(tl,i)=-cv*a_s(i,i)
         end if

!        call syncthreads()

       end do

        call syncthreads()

! Write back to global memory
      if (tx <= n .and. ty <= n) then
         a(tx,ty)=a_s(tx,ty)
      endif

     end subroutine ztrti2_gpu

! This routine computes the inverse of a complex lower non-unit triangular matrix.
! This is the blocked code that operates on a generic size matrix
 subroutine ztrtri_gpu(uplo, diag, n, a, lda, info )
      use cublas
      implicit none

      character :: uplo, diag
      integer   :: lda, n, info,sizeb
      complex(8),device:: a( lda, * )

      logical:: nounit, upper, supported
      logical, external:: lsame
      integer:: nb,nn,j,jb
      type(dim3):: threads

      if (n == 0 ) return

      upper = lsame( uplo, 'U' )
      nounit = lsame( diag, 'N' )
      IF( .not.upper .and. .not.lsame( uplo, 'L' ) ) THEN
         info = -1
      else if( .not.nounit .and. .not.lsame( diag, 'u' ) ) then
         info = -2
      else if( n.lt.0 ) theN
         info = -3
      else if( lda.lt.max( 1, n ) ) theN
         info = -5
      end iF

      ! Check if this is a supported configuratio, lower non unit triangular
      supported=(.not.upper) .and. nounit
      if ( .not.supported) then
         print *," Only L and  no-unit  implemented on GPU", upper,nounit
         info=-6
         return
      endif

      ! Block size for zttri2 on GPU is 32
       nb=32
       threads=dim3(32,32,1)

      if( nb .ge. n ) then
       ! If the problem size is smaller than block size, call directly unblocked code
           call ztrti2_gpu<<<1,threads>>>( n, a, lda )
      else
       ! Use blocked code
       !
       ! Compute inverse of lower triangular matrix
       !
            nn = ( ( n-1 ) / nb )*nb + 1
            do j = nn, 1, -nb
               jb = min( nb, n-j+1 )
               if( j+jb.le.n ) then
                !
                ! Compute rows j+jb:n of current block column
                !
                  call ztrmm( 'Left', 'Lower', 'No transpose', diag, &
                              n-j-jb+1, jb, cone, a( j+jb, j+jb ), lda, &
                              a( j+jb, j ), lda )
                  call ztrsm( 'Right', 'Lower', 'No transpose', diag, &
                              n-j-jb+1, jb, -cone, a( j, j ), lda, &
                              a( j+jb, j ), lda )
               end if
               !
               ! Compute inverse of current diagonal block
               !
               call ztrti2_gpu<<<1,threads>>>( jb, a( j, j ), lda )
            end do
      endif

     end subroutine ztrtri_gpu

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! The two following routines are the GPU implementation of the ZTRTRI LAPACK routine.
! This routine computes the inverse of a double precision lower non-unit triangular matrix.
! This is the unblocked code that operates on a block 32x32
   attributes(global) subroutine dtrti2_gpu(n,a,lda)
      use cudafor
      implicit none
      integer, value    :: n,lda
      real(8),device :: a(lda,*)
      real(8),shared :: a_s(32,32)
      real(8)        :: cv
      integer           :: tx,ty,tl,i,j,ii

     if (n == 0 ) return

      tx=threadIdx%x
      ty=threadIdx%y
      ! Linear id of the thread (tx,ty)
      tl=tx+ blockDim%x*(ty-1)

      ! Load a_d in shared memory
      if (tx <= n .and. ty <= n) then
         a_s(tx,ty)=a(tx,ty)
      endif

      call syncthreads()

     ! Compute all the diagonal elements

     if (tl <=n ) a_s(tl,tl)=one/a_s(tl,tl)

      call syncthreads()

     ! For each column working backward
      do i=n-1,1,-1

         if ( tl >i .and. tl <= n) then
            cv=zero
             do j= i+1,tl
               cv = cv + a_s(j , i ) * a_s(tl ,j)
             end do
         end if
!        call syncthreads()
         if ( tl >i .and. tl <= n) then
           a_s(tl,i)=-cv*a_s(i,i)
         end if

!        call syncthreads()

       end do

        call syncthreads()

! Write back to global memory
      if (tx <= n .and. ty <= n) then
         a(tx,ty)=a_s(tx,ty)
      endif

     end subroutine dtrti2_gpu

! This routine computes the inverse of a complex lower non-unit triangular matrix.
! This is the blocked code that operates on a generic size matrix
 subroutine dtrtri_gpu(uplo, diag, n, a, lda, info )
      use cublas
      implicit none

      character :: uplo, diag
      integer   :: lda, n, info,sizeb
      real(8),device:: a( lda, * )

      logical:: nounit, upper
      logical, external:: lsame
      integer:: nb,nn,j,jb
      type(dim3):: threads

      if (n == 0 ) return

      upper = lsame( uplo, 'U' )
      nounit = lsame( diag, 'N' )
      IF( .not.upper .and. .not.lsame( uplo, 'L' ) ) THEN
         info = -1
      else if( .not.nounit .and. .not.lsame( diag, 'u' ) ) then
         info = -2
      else if( n.lt.0 ) theN
         info = -3
      else if( lda.lt.max( 1, n ) ) theN
         info = -5
      end iF

      if (.not. upper .or. .not.nounit ) then
         print *," Only L no unit  implemented on GPU"
      endif

      ! Block size for dttri2 on GPU is 32
       nb=32
       threads=dim3(32,32,1)

      if( nb .ge. n ) then
       ! If the problem size is smaller than block size, call directly unblocked code
           call dtrti2_gpu<<<1,threads>>>( n, a, lda )
      else
       ! Use blocked code
       !
       ! Compute inverse of lower triangular matrix
       !
            nn = ( ( n-1 ) / nb )*nb + 1
            do j = nn, 1, -nb
               jb = min( nb, n-j+1 )
               if( j+jb.le.n ) then
                !
                ! Compute rows j+jb:n of current block column
                !
                  call dtrmm( 'Left', 'Lower', 'No transpose', diag, &
                              n-j-jb+1, jb, one, a( j+jb, j+jb ), lda, &
                              a( j+jb, j ), lda )
                  call dtrsm( 'Right', 'Lower', 'No transpose', diag, &
                              n-j-jb+1, jb, -one, a( j, j ), lda, &
                              a( j+jb, j ), lda )
               end if
               !
               ! Compute inverse of current diagonal block
               !
               call dtrti2_gpu<<<1,threads>>>( jb, a( j, j ), lda )
            end do
      endif

     end subroutine dtrtri_gpu

end module

#endif


text:
!##############################################################################
!# Tutorial 006f: Discretise with Q1, create a vector
!##############################################################################

module tutorial006f

  ! Include basic Feat-2 modules
  use fsystem
  use genoutput

  use triangulation
  use meshgeneration

  use element
  use spatialdiscretisation
  use linearsystemscalar

  use vectorio

  implicit none
  private

  public :: start_tutorial006f

contains

  ! ***************************************************************************

  subroutine start_tutorial006f

    ! Declare some variables.
    type(t_triangulation) :: rtriangulation
    type(t_spatialdiscretisation) :: rdiscretisation
    type(t_vectorScalar) :: rx
    character(LEN=SYS_STRLEN) :: spostdir

    ! Print a message
    call output_lbrk()
    call output_separator (OU_SEP_STAR)
    call output_line ("This is FEAT-2. Tutorial 006f")
    call output_separator (OU_SEP_MINUS)

    ! =================================
    ! Create a brick mesh
    ! =================================

    ! The mesh must always be in "standard" format.
    ! First create a 5x5-mesh on [0,1]x[0,1], then convert to standard.
    call meshgen_rectangular2DQuadMesh (rtriangulation, 0.0_DP, 1.0_DP, 0.0_DP, 1.0_DP, 4, 4)
    call tria_initStandardMeshFromRaw (rtriangulation)

    ! =================================
    ! Discretise with Q1.
    !
    ! Create a structure rdiscretisation
    ! which describes the discretisation.
    ! =================================

    call spdiscr_initDiscr_simple (rdiscretisation,EL_Q1_2D,rtriangulation)

    ! =================================
    ! Create a scalar vector
    ! =================================

    call lsyssc_createVector (rdiscretisation,rx)

    ! Fill the vector with "1.0".
    call lsyssc_clearVector (rx,1.0_DP)

    ! =================================
    ! Output to files.
    ! =================================
    call output_line ("Writing vector to text files...")

    if (sys_getenv_string("POSTDIR",spostdir)) then
      ! Write the vector to a text file.
      call vecio_writeVectorHR (rx, "vector", .true., 0, &
          trim(spostdir)//"/tutorial006f_vector.txt", "(E11.2)")

      ! Write the vector to a MATLAB file.
      call vecio_spyVector(&
          trim(spostdir)//"/tutorial006f_vector","vector",rx,.true.)

      ! Write the vector to a MAPLE file
      call vecio_writeVectorMaple (rx, "vector", .true., 0, &
          trim(spostdir)//"/tutorial006f_vector.maple", "(E11.2)")
    else
      ! Write the vector to a text file.
      call vecio_writeVectorHR (rx, "vector", .true., 0, &
          "./post/tutorial006f_vector.txt", "(E11.2)")

      ! Write the vector to a MATLAB file.
      call vecio_spyVector(&
          "./post/tutorial006f_vector","vector",rx,.true.)

      ! Write the vector to a MAPLE file
      call vecio_writeVectorMaple (rx, "vector", .true., 0, &
          "./post/tutorial006f_vector.maple", "(E11.2)")
    end if

    ! =================================
    ! Cleanup
    ! =================================

    ! Release the vector
    call lsyssc_releaseVector (rx)

    ! Release the Q1-discretisation
    call spdiscr_releaseDiscr (rdiscretisation)

    ! Release the triangulation
    call tria_done (rtriangulation)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense