Numerical Methods 資工所碩一 r99922026 王裕博
1. My code:
#include
#include
#include
#include
#include
#include
#include
static sigjmp_buf jmpbuf;
void handler(int sig, siginfo_t *sip,void *uap){
unsigned code,addr;
code = sip->si_code;
addr = (unsigned) sip->si_addr;
switch (code){
case FPE_FLTOVF: fprintf(stderr , "overflow\n"); break;
case FPE_FLTUND: fprintf(stderr , "underflow\n"); break;
case FPE_FLTDIV: fprintf(stderr,"divided by zero\n"); break;
case FPE_FLTRES: fprintf(stderr,"inexact\n"); break;
case FPE_FLTINV: fprintf(stderr,"invaild operation\n"); break;
}
fprintf(stderr , "fp floating point exception %x at address %x \n", code , addr);
feenableexcept(FE_OVERFLOW| FE_DIVBYZERO| FE_INVALID);
siglongjmp(jmpbuf,1);
}
int main() {
double x;
struct sigaction act,oldact;
act.sa_handler = handler;
act.sa_flags = SA_SIGINFO;
if( sigaction(SIGFPE,&act,&oldact) != 0 ){
fprintf(stderr, "sigaction error\n");
return 0 ;
}
feclearexcept(FE_ALL_EXCEPT);
feenableexcept(FE_OVERFLOW| FE_DIVBYZERO| FE_INVALID);
x = DBL_MIN;
printf("min_normal = %g\n", x);
if(sigsetjmp(jmpbuf,1)==0) x/= 13.0;
printf("min_normal / 13.0 = %g\n", x);
x = DBL_MAX;
printf("max_normal = %g\n", x);
if(sigsetjmp(jmpbuf,1)==0){
x = x*x;
}
printf("max_normal * max_normal = %g\n",x);
return 0;
}
In this simulation, I used “feclearexcept()” function to clear all signal, and used
“feenableexcept()” function to open the “common” signal(overflow, divided by zero,
invalid operation ). And I used “sigaction()” to replace the function which be used
when signal happened with my own function “handler”. In function handler, I just do
what the course slide showed. Finally, the “sigsetjmp()” and “siglongjmp()” deal with
the infinite production of the signal. And results are as bellows:
r99922026@linux14 [~/nm] gcc -lm test.c
test.c: In function 'handler':
test.c:13:8: warning: cast from pointer to integer of different size
test.c: In function 'main':
test.c:28:16: warning: assignment from incompatible pointer type
r99922026@linux14 [~/nm] ./a.out
min_normal = 2.22507e-308
min_normal / 13.0 = 1.7116e-309
max_normal = 1.79769e+308
overflow
fp floating point exception 4 at address 400a0f
max_normal * max_normal = 1.79769e+308
which is similar to the course slides.
When underflow happened, because we don’t set the flag, it will go to the default
road, and just treated it as a minimum number. But things changed when overflow
happened. It will go to our handler and print the message what we want to print!
2. Here is my code
cholevel1.m
function d = cohlevel1(c)
d = c;
for k = 1:2500
d(k,k) = sqrt(d(k,k));
d(k+1:2500,k) = d(k+1:2500,k) / d(k,k);
d(k+1:2500,k+1:2500) = d(k+1:2500,k+1:2500)-(d(k+1:2500,k) *
d(k+1:2500,k)');
end
cholevel2.m
function d = cohlevel2(c)
d = c;
for k = 1:2500
d(k,k) = sqrt(d(k,k));
d(k+1:2500,k) = d(k+1:2500,k) / d(k,k);
for j = k+1:2500
d(j:2500,j) = d(j:2500,j)-(d(j:2500,k)*d(j,k));
end
end
cholevel3.m
function d = cohlevel3(c)
d = c;
for k = 1:2500
d(k,k) = sqrt(d(k,k));
d(k+1:2500,k) = d(k+1:2500,k) / d(k,k);
for j = k+1:2500
for i=j:2500
d(i,j) = d(i,j) - (d(i,k)*d(j,k));
end
end
end
and I run it this way:
>> clear all
>> a = rand(2500,2500);
>> b = a';
>> c = a*b;
>> t0 = cputime;d = cholevel1(c);cpu_time = cputime-t0
cpu_time =
205.4221
>> t0 = cputime;d = cholevel2(c);cpu_time = cputime-t0
cpu_time =
53.4459
>> t0 = cputime;d = cholevel3(c);cpu_time = cputime-t0
cpu_time =
86.1750
>> t0 = cputime;d = chol(c);cpu_time = cputime-t0
cpu_time =
1.1076
And we can see that:
Level 1 time > Level 3 time > Level 2 time > chol() time.