fortran - Prevent changing variables with intent(in) -
so reading following question (correct use of fortran intent() large arrays) learned defining variable intent(in) isn't enough, since when variable passed subroutine/function, can changed again. how can avoid this? in original thread talked putting subroutine module, doesn't me. example want calculate determinant of matrix lu-factorization. therefore use lapack function zgetrf, function alters input matrix , compiler don't displays warnings. can do?
module mathelper implicit none contains subroutine initmat(aa) real*8 :: u double complex, dimension(:,:), intent(inout) :: aa integer :: row, col, counter counter = 1 row=1,size(aa,1) col=1,size(aa,2) aa(row,col)=cmplx(counter ,0) counter=counter+1 end end end subroutine initmat !subroutine write matrix file !input: aa - double complex matrix ! fid - integer file id ! fname - file name ! stat - integer status =replace[0] or old[1] subroutine writemat(aa,fid, fname, stat) integer :: fid, stat character(len=*) :: fname double complex, dimension(:,:), intent(in) :: aa integer :: row, col character (len=64) :: fmtstring !opening file given options if(fid /= 0) if(stat == 0) open(unit=fid, file=fname, status='replace', & action='write') else if(stat ==1) open(unit=fid, file=fname, status='old', & action='write') else print*, 'error while trying open file id', fid return end if end if !initializing matrix print format write(fmtstring,'(i0)') size(aa,2) fmtstring = '('// trim(fmtstring) //'("{",es10.3, ",", 1x, es10.3,"}",:,1x))' !write(*,*) fmtstring !writing matrix file iterating through each row row=1,size(aa,1) write(fid,fmt = fmtstring) aa(row,:) enddo write(fid,*) '' end subroutine writemat !function calculate determinant of input !input: aa - double complex matrix !output determinantmat - double complex, ! 0 if aa not square matrix function determinantmat(aa) double complex, dimension(:,:), intent(in) :: aa double complex :: determinantmat integer, dimension(min(size(aa,1),size(aa,2)))& :: ipiv integer :: ii, info !check if not square matrix, set determinant 0 if(size(aa,1)/= size(aa,2)) determinantmat = 0 return end if !compute lu facotirzation lapack function call zgetrf(size(aa,1),size(aa,2), aa,size(aa,1), ipiv,info) if(info /= 0) determinantmat = cmplx(0.d0, 0.d0) return end if determinantmat = cmplx(1.d0, 0.d0) !determinant of triangular matrix product of diagonal elements ii=1,size(aa,1) if(ipiv(ii) /= ii) !a permutation done, factor of -1 determinantmat = -determinantmat *aa(ii,ii) else !no permutation, no -1 determinantmat = determinantmat*aa(ii,ii) end if end end function determinantmat end module mathelper !*********************************************************************** !module stores matrix elements, dimension, trace, determinant program test use mathelper implicit none double complex, dimension(:,:), allocatable :: aa, bb integer :: n, fid fid = 0; allocate(aa(3,3)) call initmat(aa) call writemat(aa,0,' ', 0) print*, 'determinante: ',determinantmat(aa) !changes aa call writemat(aa,0, ' ', 0) end program test ps: using ifort compiler v15.0.3 20150407
i not have ifort @ home, may want try compiling '-check interfaces' , maybe '-ipo'. may need path 'zgetrf' '-check interfaces' work, , if not source may not help. if declare 'function determinantmat' 'pure function determinantmat' pretty sure complain because 'zgetrf' not known pure nor elemental. try ^this stuff^ first.
if lapack has module, zgetrf known be, or not be, pure/elemental. https://software.intel.com/en-us/articles/blas-and-lapack-fortran95-mod-files
i suggest add compile line:
-check interfaces -ipo during initial build (take out speed once works):
-check -warn making temporary array 1 way around it. (i have not compiled this, conceptual exemplar.)
pure function determinantmat(aa) use lapack95 !--new line--! implicit none !--new line--! double complex, dimension(:,:) , intent(in ) :: aa double complex :: determinantmat !<- output !--internals-- integer, dimension(min(size(aa,1),size(aa,2))) :: ipiv !!--next line new-- double complex, dimension(size(aa,1),size(aa,2)) :: aa_temp !!<- have no idea if work, may need allocatable?? integer :: ii, info !check if not square matrix, set determinant 0 if(size(aa,1)/= size(aa,2)) determinantmat = 0 return end if !compute lu factorization lapack function !!--next line new-- aa_temp = aa !--initialise aa_temp same aa--! call zgetrf(size(aa_temp,1),size(aa_temp,2), aa_temp,size(aa_temp,1), ipiv,info) if(info /= 0) determinantmat = cmplx(0.d0, 0.d0) return end if determinantmat = cmplx(1.d0, 0.d0) !determinant of triangular matrix product of diagonal elements ii=1,size(aa_temp,1) if(ipiv(ii) /= ii) !a permutation done, factor of -1 determinantmat = -determinantmat *aa_temp(ii,ii) else !no permutation, no -1 determinantmat = determinantmat*aa_temp(ii,ii) end if end end function determinantmat with 'use lapack95' not need pure, if wanted pure want explicitly so.
Comments
Post a Comment